Solar Magnetohydrodynamics: Solving Magnetic Diffusion in a Sphere

The Sun’s Magnetic Field as a PDE Problem

The sun is a magnetised plasma — a turbulent, rotating, electrically conducting fluid threaded by magnetic fields of extraordinary complexity. Sunspots, solar flares, coronal mass ejections, and the 11-year activity cycle all trace back to the evolution of this field. At the mathematical core of it all sits a deceptively clean PDE: the **MHD induction equation**.

In this post I solve two instances of this equation in spherical geometry using **Dedalus**, a Python spectral PDE library. The first simulation evolves the magnetic field on the solar surface — a 2D problem on the sphere $S^2$. The second evolves the field through the solar interior — a 3D problem inside a conducting ball. Both produce animated GIFs that reveal the underlying physics.

The MHD Induction Equation

In a conducting fluid, the magnetic field $\mathbf{B}$ evolves according to the **induction equation**:

$$\frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{u} \times \mathbf{B}) + \eta \, \nabla^2 \mathbf{B}$$

$(1)$

The first term on the right is the **inductive term**: the velocity field $\mathbf{u}$ stretches, twists, and amplifies the magnetic field. The second term is **Ohmic diffusion**: resistivity $\eta$ causes the field to smooth out and decay. The ratio of these two effects is the **magnetic Reynolds number** $R_m = UL/\eta$, which controls whether the field is amplified (large $R_m$, as in the solar dynamo) or simply diffuses away (small $R_m$).

When the flow vanishes ($\mathbf{u} = 0$), equation (1) reduces to **pure magnetic diffusion**:

$$\frac{\partial \mathbf{B}}{\partial t} = \eta \, \nabla^2 \mathbf{B}$$

$(2)$

This is vector heat equation in $\mathbf{B}$. The field decays through its **free decay modes** — eigenfunctions of the Laplacian in the given geometry, each decaying exponentially at a rate set by $\eta$ and the spatial scale. In a sphere, these modes are expressed in terms of spherical harmonics, which is precisely where spectral methods become natural.

For our simulations we model $B_r$ — the radial (normal) component of the magnetic field — as a scalar proxy. This is the component that most directly controls the force the field exerts on the overlying atmosphere and is the quantity mapped by solar observatories.

Spectral Methods on a Sphere

On a periodic rectangular domain, Fourier modes are the natural basis. On a sphere the analogous basis is **spherical harmonics** $Y_\ell^m(\theta, \phi)$, which are eigenfunctions of the Laplace–Beltrami operator:

$$\nabla^2_{S^2} Y_\ell^m = -\frac{\ell(\ell+1)}{R^2} Y_\ell^m$$

$(3)$

Each harmonic degree $\ell$ corresponds to a spatial scale $\sim \pi R / \ell$. High-$\ell$ modes are small-scale structure; $\ell = 1$ is the global dipole. The diffusion equation in spectral space becomes independent equations for each mode, each decaying as $e^{-\eta \ell(\ell+1) t / R^2}$ — high-$\ell$ modes die first, leaving only the large-scale dipole at late times.

Dedalus 3 implements this in two forms:

  • `SphereBasis` — for PDEs on the 2D sphere surface $S^2$ (colatitude + azimuth)
  • `BallBasis` — for PDEs in the 3D ball $r \leq R$ (radius + colatitude + azimuth)

Both use spin-weighted spherical harmonic transforms internally, which correctly handle the coordinate singularities at the poles and at $r = 0$.

Simulation 1 — Surface Flux Transport on $S^2$

The **surface flux transport** (SFT) model is a workhorse of solar physics. It evolves $B_r$ on the photosphere under the combined effects of differential rotation, meridional circulation, and supergranular diffusion. Here we isolate the diffusion component, which is dominant over long time scales:

$$\frac{\partial B_r}{\partial t} = \eta \, \nabla^2_{S^2} B_r$$

$(4)$

The initial condition mimics a realistic solar disc: three **bipolar active regions** (sunspot groups) in each hemisphere, with the correct Hale polarity law — northern and southern hemisphere pairs have opposite leading polarities. The regions are Gaussian blobs of opposite sign, placed at active latitudes ($\pm 20$–$40°$).

Python 3.13 — Surface Flux Transport: Magnetic Diffusion on S²

import os
os.environ["OMP_NUM_THREADS"] = "1"

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import dedalus.public as d3
import logging

logging.getLogger('dedalus').setLevel(logging.ERROR)

# Parameters
Nphi, Ntheta = 128, 64
eta          = 0.015        # surface magnetic diffusivity
stop_time    = 2.5
timestep     = 0.002        # fine step for accuracy

# 2D sphere surface basis (spin-weighted spherical harmonics)
coords = d3.S2Coordinates('phi', 'theta')
dist   = d3.Distributor(coords, dtype=np.float64)
basis  = d3.SphereBasis(coords, (Nphi, Ntheta), np.float64,
                        radius=1, dealias=(3/2, 3/2))

B = dist.Field(name='B', bases=basis)

# PDE: surface diffusion — Laplace-Beltrami operator on S²
problem = d3.IVP([B], namespace=locals())
problem.add_equation("dt(B) = eta * lap(B)")

solver = problem.build_solver(d3.RK443)
solver.stop_sim_time = stop_time

# Initial condition: bipolar active regions with Hale polarity law
phi_g, theta_g = dist.local_grids(basis)
lat = np.pi/2 - theta_g   # colatitude → latitude

def bipole(lat0, lon0, size=25, sep=0.35, amp=1.0):
    pos = amp  * np.exp(-size * ((lat - lat0)**2 + (phi_g - lon0)**2))
    neg = -amp * np.exp(-size * ((lat - lat0)**2 + (phi_g - lon0 - sep)**2))
    return pos + neg

B['g'] = (
    bipole( 0.35, 1.0, size=30, sep= 0.30, amp=1.0) +   # NH pair 1
    bipole( 0.25, 3.5, size=25, sep= 0.40, amp=0.8) +   # NH pair 2
    bipole( 0.40, 5.5, size=20, sep= 0.35, amp=0.9) +   # NH pair 3
    bipole(-0.35, 1.2, size=30, sep=-0.30, amp=-1.0) +  # SH pair 1 (reversed)
    bipole(-0.25, 4.0, size=25, sep=-0.40, amp=-0.8) +  # SH pair 2
    bipole(-0.40, 5.8, size=20, sep=-0.35, amp=-0.9)    # SH pair 3
)

# Non-linear frame sampling — dense early where structure evolves fastest
frame_times = np.concatenate([
    np.linspace(0,   0.8, 70),   # 70 frames in first 0.8 time units
    np.linspace(0.8, 2.5, 30),   # 30 frames for slow tail
])
frames = []
frame_idx = 0

while solver.proceed and frame_idx < len(frame_times):
    if solver.sim_time >= frame_times[frame_idx]:
        frames.append(B['g'].copy())
        frame_idx += 1
    solver.step(timestep)
Figure 1
Figure 1. Solar surface magnetic field evolving under diffusion on the 2D sphere $S^2$. Red: positive $B_r$ (north polarity). Blue: negative $B_r$ (south polarity). Dashed line: solar equator. The bipolar active region pairs spread and eventually cancel across the equator, leaving a weak polar field — the same process that reverses the sun’s polar fields between cycles.

The simulation reveals the key physics: active regions with opposite polarities spread toward the equator, where they meet and cancel. The higher-latitude flux preferentially migrates toward the poles. Over a full solar cycle timescale, this process reverses the sun’s polar magnetic field — the mechanism underlying Hale’s polarity law.

Simulation 2 — Interior Magnetic Diffusion in a 3D Ball

Inside the sun, the magnetic field diffuses through a conducting medium according to equation (2). We solve this in 3D using `BallBasis`, which covers the full interior of a sphere from the centre $r = 0$ to the outer radius $r = R$.

The boundary condition at $r = R$ is chosen as $B_r = 0$, modelling a perfectly conducting outer shell (the lowest-order approximation to matching onto an external potential field). Regularity at $r = 0$ is enforced automatically by the basis.

The initial condition is a **dipolar field** — the $\ell = 1$, $m = 0$ structure that dominates the sun’s large-scale field:

$$T(r, \theta, \phi, 0) = r(1 – r^2)\cos\theta + \epsilon \, r(1-r^2)\sin^2\theta \cos 2\phi$$

$(5)$

The small $\epsilon$ perturbation breaks axial symmetry, making the diffusion more visually interesting.

Python 3.13 — Interior Magnetic Diffusion in a 3D Ball (Meridional Cross-Section)

import os
os.environ["OMP_NUM_THREADS"] = "1"

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import dedalus.public as d3
import logging

logging.getLogger('dedalus').setLevel(logging.ERROR)

# Parameters
Nphi, Ntheta, Nr = 64, 32, 32
R        = 1.0
eta      = 0.08
stop_time = 3.0
timestep  = 0.004        # fine step for accuracy

# 3D ball basis (full sphere interior)
coords = d3.SphericalCoordinates('phi', 'theta', 'r')
dist   = d3.Distributor(coords, dtype=np.float64)
ball   = d3.BallBasis(coords, (Nphi, Ntheta, Nr), np.float64,
                      radius=R, dealias=(3/2, 3/2, 3/2))
S2     = ball.surface   # outer 2D sphere for boundary condition

# Fields + tau for boundary condition enforcement
T   = dist.Field(name='T',   bases=ball)
tau = dist.Field(name='tau', bases=S2)

lift_basis = ball.derivative_basis(1)
lift       = lambda A: d3.Lift(A, lift_basis, -1)

# PDE: 3D magnetic diffusion with Dirichlet BC at r=R
problem = d3.IVP([T, tau], namespace=locals())
problem.add_equation("dt(T) - eta*lap(T) + lift(tau) = 0")
problem.add_equation("T(r=1) = 0")   # perfectly conducting outer shell

solver = problem.build_solver(d3.RK443)
solver.stop_sim_time = stop_time

# Initial condition: dipole + small non-axisymmetric perturbation
phi_g, theta_g, r_g = dist.local_grids(ball)
T.change_scales(1)
T['g'] = (r_g * (1 - r_g**2) * np.cos(theta_g)
          + 0.3 * r_g * (1 - r_g**2) * np.sin(theta_g)**2 * np.cos(2*phi_g))

# Non-linear frame sampling — dense early where structure changes rapidly
frame_times = np.concatenate([
    np.linspace(0,   1.0, 65),   # 65 frames in first 1.0 time units
    np.linspace(1.0, 3.0, 25),   # 25 frames for slow decay tail
])
frames = []
frame_idx = 0

while solver.proceed and frame_idx < len(frame_times):
    if solver.sim_time >= frame_times[frame_idx]:
        T.change_scales(1)
        frames.append(T['g'].mean(axis=0).copy())   # (Ntheta, Nr)
        frame_idx += 1
    solver.step(timestep)

The animation shows the meridional cross-section of the sphere (the equatorial plane runs horizontally, the polar axis vertically). Red is positive field, blue is negative.

Figure 2
Figure 2. Interior magnetic diffusion in a 3D conducting sphere, shown as a meridional cross-section. The initial dipolar structure — positive (red) in the upper hemisphere, negative (blue) in the lower — decays through Ohmic diffusion. Small-scale structure dies first; the global dipole ($\ell=1$ mode) persists longest, exactly as predicted by the spherical harmonic decay theory.

The simulation confirms the theoretical prediction: higher-$\ell$ modes (finer spatial structure) decay faster, leaving the smooth dipolar mode as the long-lived remnant. This is the free decay of the geomagnetic field — relevant to understanding why the Earth’s dipole weakens and occasionally reverses on geological timescales.

Figure 3 — 3D Rotating Sphere Magnetic Diffusion
Figure 3. Magnetic field intensity on the sphere surface at r ≈ 0.76, rotating 540° — jet colormap, same simulation as Figure 2.

A note on notation: Figures 1 and 2 label their colorbars Br [a.u.] — the radial component of the magnetic field, the quantity that points outward or inward through the sphere surface and is directly measured in solar magnetograms. Figure 3 labels its colorbar T [a.u.] because the 3D ball solver uses a scalar field T as a proxy. The governing equation is identical — ∂T/∂t = η∇²T — so T and Br obey exactly the same diffusion physics. The different name simply reflects that the spectral solver treats the unknown as a generic scalar; interpreting it as the radial magnetic field component is the physical assignment we impose through the initial condition and boundary conditions.

Where These Equations Appear

**Solar physics.** The surface flux transport model is used operationally at observatories worldwide to forecast how active region flux will redistribute across the photosphere and ultimately drive the solar cycle. The 3D induction equation governs the solar dynamo — the engine that amplifies and sustains the field against Ohmic decay.

**Space weather.** The rate at which magnetic flux reaches the polar caps determines the intensity of geomagnetic storms. Accurate SFT modelling directly improves 27-day solar wind forecasts.

**Planetary interiors.** The Earth’s magnetic field is generated in the liquid outer core. The same free-decay modes govern how long a frozen-in fossil field would survive — a key question in paleomagnetism and planetary science.

**Tokamak plasma confinement.** In fusion reactors, magnetic field topology must be controlled precisely. Resistive MHD instabilities — tearing modes — are governed by variants of the same diffusion equation in toroidal geometry.

**Neutron stars.** A newly formed neutron star has an enormous magnetic field ($\sim 10^{12}$ G). As the crust solidifies, Ohmic decay determines how rapidly this field evolves over millions of years — the same physics, at fantastically different scales.

Takeaway

Spectral methods on spherical domains are not just elegant — they are the natural choice for problems where the geometry is spherical and the solutions are smooth. Dedalus makes it possible to implement these solvers in Python with minimal boilerplate: the PDE is written symbolically, the library handles the transforms, and the results are exact to spectral precision.

The solar magnetic field is a system where this matters. One extra harmonic degree of resolution can mean the difference between capturing a key instability and missing it entirely.

Need a reliable numerical implementation for your problem? Get in touch.