Visualizing the Earth's dipolar magnetic field

The magnetic field due to a magnetic dipole moment, $\boldsymbol{m}$ at a point $\boldsymbol{r}$ relative to it may be written

\boldsymbol{B}(\boldsymbol{r}) = \frac{\mu_0}{4\pi r^3}[3\boldsymbol{\hat{r}(\boldsymbol{\hat{r}} \cdot \boldsymbol{m}) - \boldsymbol{m}}],

where $\mu_0$ is the vacuum permeability. In geomagnetism, it is usual to write the radial and angular components of $\boldsymbol{B}$ as:

\begin{align*} B_r & = -2B_0\left(\frac{R_\mathrm{E}}{r}\right)^3\cos\theta, \\ B_\theta & = -B_0\left(\frac{R_\mathrm{E}}{r}\right)^3\sin\theta, \\ B_\phi &= 0, \end{align*}

where $\theta$ is polar (colatitude) angle (relative to the magnetic North pole), $\phi$ is the azimuthal angle (longitude), and $R_\mathrm{E}$ is the Earth's radius, about 6370 km. See below for a derivation of these formulae.

With this definition, $$B_0$$ denotes the magnitude of the mean value of the field at the magnetic equator on the Earth's surface, about $31.2\; \mathrm{\mu T}$ . With these definitions, we can plot the dipole component of the Earth's magnetic field as a Matplotlib streamplot. We first construct a meshgrid of $$(x,y)$$ coordinates and convert them into polar coordinates with the relations:

\begin{align*} r &= \sqrt{x^2 + y^2},\\ \theta &= \mathrm{atan2}(y/x), \end{align*}

where the two-argument arctangent function, $\mathrm{atan2}$ , is implemented in NumPy as arctan2. We can then use the formulae above to calculate $$B_r$$ and $B_\theta$ and convert back to Cartesian coordinates to plot the streamplot. The $\theta$ coordinate is offset by $\alpha = 9.6^\circ$ to account for the current tilt of the Earth's magnetic dipole with respect to its rotational axis.

Earth's dipolar magnetic field

import sys
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle

# Mean magnitude of the Earth's magnetic field at the equator in T
B0 = 3.12e-5
# Radius of Earth, Mm (10^6 m: mega-metres!)
RE = 6.370
# Deviation of magnetic pole from axis
alpha = np.radians(9.6)

def B(r, theta):
    """Return the magnetic field vector at (r, theta)."""
    fac = B0 * (RE / r)**3
    return -2 * fac * np.cos(theta + alpha), -fac * np.sin(theta + alpha)

# Grid of x, y points on a Cartesian grid
nx, ny = 64, 64
XMAX, YMAX = 40, 40
x = np.linspace(-XMAX, XMAX, nx)
y = np.linspace(-YMAX, YMAX, ny)
X, Y = np.meshgrid(x, y)
r, theta = np.hypot(X, Y), np.arctan2(Y, X)

# Magnetic field vector, B = (Ex, Ey), as separate components
Br, Btheta = B(r, theta)
# Transform to Cartesian coordinates: NB make North point up, not to the right.
c, s = np.cos(np.pi/2 + theta), np.sin(np.pi/2 + theta)
Bx = -Btheta * s + Br * c
By = Btheta * c + Br * s

fig, ax = plt.subplots()

# Plot the streamlines with an appropriate colormap and arrow style
color = 2 * np.log(np.hypot(Bx, By))
ax.streamplot(x, y, Bx, By, color=color, linewidth=1, cmap=plt.cm.inferno,
              density=2, arrowstyle='->', arrowsize=1.5)

# Add a filled circle for the Earth; make sure it's on top of the streamlines.
ax.add_patch(Circle((0,0), RE, color='b', zorder=100))

ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_xlim(-XMAX, XMAX)
ax.set_ylim(-YMAX, YMAX)
ax.set_aspect('equal')
plt.show()

Derivation

For now, consider the dipole to lie along the $$y$$ -axis: $\boldsymbol{m} = m\boldsymbol{\hat{\jmath}}$ and define $\boldsymbol{B_0}$ to be the field perpendicular to $\boldsymbol{m}$ at a distance $R_\mathrm{E}$ from it (i.e. on the Earth's magnetic equator):

\boldsymbol{B_0} = -\frac{\mu_0m}{4\pi R_\mathrm{E}^3} \boldsymbol{\hat{\jmath}} = B_0\boldsymbol{\hat{\jmath}}.

Then,

\boldsymbol{B}(\boldsymbol{r}) = -B_0 \left(\frac{R_\mathrm{E}}{r}\right)^3 [3\cos\theta \boldsymbol{\hat{r}} - \boldsymbol{\hat{\jmath}}].

The radial component of the magnetic field is therefore

B_r = \boldsymbol{B}\cdot\boldsymbol{\hat{r}} = -B_0\left(\frac{R_\mathrm{E}}{r}\right)^3[3\cos\theta - \cos\theta] = -2B_0\left(\frac{R_\mathrm{E}}{r}\right)^3\cos\theta,

Its angular component is perpendicular to this:

B_\theta = \boldsymbol{B}\cdot\boldsymbol{\hat{\theta}} = - B_0 \left(\frac{R_\mathrm{E}}{r}\right)^3 [-\boldsymbol{\hat{\jmath}} \cdot \boldsymbol{\hat{\theta}}] = -B_0\left(\frac{R_\mathrm{E}}{r}\right)^3\sin\theta.

Finally, $B_\phi = 0$ : the field is symmetric about the axis of the dipole.

Visualizing the Earth's dipolar magnetic field

Derivation

Further Reading