Plasma and MHD — Magnetohydrodynamics

1. What is Plasma: Ionised Gas

Plasma forms when a gas is heated so intensely that electrons are stripped from atoms, creating a mixture of free electrons and ions. Key parameters:

The Debye length is the distance over which electric fields are screened. Beyond this radius the plasma is quasi-neutral: n_e ≈ n_i. In the solar corona λ_D ~ 10⁻² m, in laboratory plasma ~ 10⁻⁴ m.

2. Lorentz Force and Charged Particle Motion

A charge q in an electric and magnetic field experiences the Lorentz force:

Larmor Gyration

In a uniform magnetic field B = B_zẑ, a particle moves in a circle (in the xy plane) and rectilinearly along z. Cyclotron (Larmor) frequency and radius:

E×B Drift

In the presence of both fields the particle drifts perpendicular to both E and B:

3. MHD Equations: Navier-Stokes + Maxwell

MHD describes plasma as an ideal or resistive conducting fluid. Single-fluid MHD is a simplified model where electrons and ions are not distinguished separately. The equations are:

For ideal MHD η = 0: the field is frozen into the fluid ("frozen-in flux theorem"). Magnetic field lines move together with the plasma — a fundamental property.

Pressure / Energy Equation

The parameter β = p/p_B = 2μ₀p/B² defines the regime: β ≪ 1 — magnetically dominated (corona, thin disc); β ~ 1 — MHD; β ≫ 1 — hydrodynamic.

4. Pinch Effect and Magnetic Confinement

The pinch (Z-pinch) is perhaps the simplest MHD configuration. A current J_z flows through the plasma, generating an azimuthal field B_φ. The Ampère force J × B is directed toward the axis — the plasma is compressed. Equilibrium:

Pinch Instabilities

5. Rayleigh-Taylor Instability in Plasma

The Rayleigh-Taylor (RT) instability arises when a heavier fluid lies above a lighter one in a gravitational field — or in MHD when dense plasma is held above a rarefied region by a magnetic field.

Dispersion Relation

The magnetic field stabilises short-wavelength perturbations, but long-wavelength ones (k < k_c) remain unstable. In laser-driven fusion, the RT instability is the main obstacle to uniform target compression.

Magnetic Reconnection

When two antiparallel field configurations approach each other, a thin layer where B ≈ 0 (a current sheet) can tear — the field reorganises ("reconnects"), releasing enormous kinetic and thermal energy. Solar flares are exactly this phenomenon at a scale of ~10⁹ m.

6. Alfvén Waves

An Alfvén wave is the MHD analogue of a wave on a taut string. The magnetic tension force B²/μ₀ acts as the "spring", ρ as the "mass per unit length". Wave speed:

Types of MHD Waves

Alfvén waves heat the solar corona, transport energy from the convection zone to the corona, and play a role in accelerating the solar wind.

7. Solar Wind and Earth's Magnetosphere

Parker (1958) showed that the hot solar corona cannot be in equilibrium with interstellar space — it must expand. The solution of MHD equations through the critical point r_c = GM_⊙/(2c_s²):

Interaction with the Magnetosphere

The solar wind encounters the geomagnetic field and forms: a bow shock at ~14 R_⊕ where the wind is decelerated; the magnetopause — the boundary where wind pressure equals magnetic field pressure; the magnetotail — the nightside tail extending ~100 R_⊕.

8. Simple MHD Simulation: JS Implementation

We implement 2D ideal MHD on an explicit FD scheme (Lax-Friedrichs). System state: (ρ, v_x, v_y, B_x, B_y, p) at each grid node.

class MHD2D {
  constructor(Nx, Ny, dx = 1, dy = 1) {
    this.Nx = Nx; this.Ny = Ny;
    this.dx = dx; this.dy = dy;
    this.gamma = 5 / 3;
    const N = Nx * Ny;
    this.rho = new Float64Array(N).fill(1); // density
    this.vx  = new Float64Array(N);
    this.vy  = new Float64Array(N);
    this.Bx  = new Float64Array(N);
    this.By  = new Float64Array(N);
    this.p   = new Float64Array(N).fill(1); // pressure
  }

  alfvenSpeed(i, j) {
    const k = i * this.Ny + j;
    const B2 = this.Bx[k]**2 + this.By[k]**2;
    return Math.sqrt(B2 / this.rho[k]); // μ₀ = 1 in code units
  }

  cflTimestep() {
    let maxSpeed = 0;
    for (let i = 0; i < this.Nx; i++) {
      for (let j = 0; j < this.Ny; j++) {
        const k = i * this.Ny + j;
        const cs = Math.sqrt(this.gamma * this.p[k] / this.rho[k]);
        const va = this.alfvenSpeed(i, j);
        const v  = Math.hypot(this.vx[k], this.vy[k]);
        maxSpeed = Math.max(maxSpeed, v + cs + va);
      }
    }
    return 0.4 * Math.min(this.dx, this.dy) / maxSpeed; // CFL = 0.4
  }

  initOrzagTang() {
    // Classic test: Orszag-Tang vortex
    const { Nx, Ny } = this;
    for (let i = 0; i < Nx; i++) {
      for (let j = 0; j < Ny; j++) {
        const k = i * Ny + j;
        const x = 2 * Math.PI * i / Nx;
        const y = 2 * Math.PI * j / Ny;
        this.rho[k] = 25 / (36 * Math.PI);
        this.vx[k] = -Math.sin(y);
        this.vy[k]  = Math.sin(x);
        this.Bx[k] = -Math.sin(y) / (4 * Math.PI);
        this.By[k]  = Math.sin(2 * x) / (4 * Math.PI);
        this.p[k]  = 5 / (12 * Math.PI);
      }
    }
  }
}

Visualisation: Density Colour Map + B Field Arrows

function renderMHD(canvas, mhd) {
  const ctx = canvas.getContext('2d');
  const img = ctx.createImageData(mhd.Nx, mhd.Ny);
  const maxRho = Math.max(...mhd.rho);
  const minRho = Math.min(...mhd.rho);

  for (let i = 0; i < mhd.Nx; i++) {
    for (let j = 0; j < mhd.Ny; j++) {
      const k = i * mhd.Ny + j;
      const t = (mhd.rho[k] - minRho) / (maxRho - minRho);
      const px = (j * mhd.Nx + i) * 4;
      img.data[px]   = lerp(10, 255, t);   // R: blue→hot
      img.data[px+1] = lerp(10, 80, t);
      img.data[px+2] = lerp(180, 10, t);
      img.data[px+3] = 255;
    }
  }
  ctx.putImageData(img, 0, 0);
}