The Hasegawa-Wakatani drift-wave model

Adam Zettel · September 25, 2026 · written with an agent

ω = \frac{κ k _{y}}{1 + k _{⊥}^{2}}, \frac{d}{d t} (ζ n) = (- C / k^{2} - C / k^{2} + iκ k_{y} / k^{2} - C - C) (ζ n)

The Hasegawa-Wakatani equations are the simplest model of drift-wave turbulence at the edge of a tokamak that contains the instability driving it. Five rungs build their linear theory: the drift-wave frequency, the growth rate as the eigenvalue of a 2×2 matrix (0.098684 at $C = 1$ , as numpy's eigensolver gives), the return to the drift wave as $C \to \infty$ at the predicted rates $1/ C$ and $1/ C^{2}$ , a zonal flow that keeps $k^{2} / (1 + k^{2})$ of its start as it relaxes, and the one-line change in the modified model that leaves it untouched.

Open in the workbench Download .vivOutputs below are from the site build’s run; Run all runs it again in your browser.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

The edge of a magnetically confined plasma has a steep density gradient, and the gradient drives drift waves: ripples in density and electric potential that travel across the magnetic field, along the surfaces of constant density. The Hasegawa-Wakatani equations (1983) describe them with two fields in the plane perpendicular to $B$ : the density fluctuation $n$ and the electric potential $ϕ$ , whose Laplacian $ζ = \nabla^{2} ϕ$ is the vorticity of the $E \times B$ flow.

\partial_{t} ζ + [ϕ, ζ] = C (ϕ - n), \partial_{t} n + [ϕ, n] = - κ \partial_{y} ϕ + C (ϕ - n) .

$κ$ is the density gradient, the drive. $C$ is the adiabaticity: how fast electrons stream along the field line to keep the density in step with the potential. Everything is dimensionless, with lengths in units of the ion sound gyroradius $ρ_{s}$ , and here $κ = 1$ . The brackets $[ϕ, \cdot]$ are the nonlinear advection; every rung here is linear, so they drop out.

posts/refs/hasegawa_wakatani.py computes the reference values with numpy.linalg.eigvals and scipy.linalg.expm.

The drift wave

In the limit $C \to \infty$ the electrons follow the potential at once and $n = ϕ$ . Subtracting the two equations then leaves one, the Hasegawa-Mima equation $\partial_{t} (\nabla^{2} ϕ - ϕ) = κ \partial_{y} ϕ$ (linear part). A wave $e^{i (k_{x} x + k_{y} y - ω t)}$ solves it when

$ω = \frac{κ k _{y}}{1 + k _{x}^{2} + k _{y}^{2}} .$

The wave moves only in $y$ , along the surfaces of constant density. The 1 in the denominator comes from $n = ϕ$ ; the $k^{2}$ comes from the ion polarization drift, which is stronger for small eddies. So the frequency first rises with $k_{y}$ , peaks at $k_{y} = 1$ (for $k_{x} = 0$ ), and then falls.

new: omega_hm

Rung 1 tests 3 of 3 pass

the drift wave at k = (0, 1) has frequency 1/2 (passes)
at k = (0.3, 0.7) it is 0.7 / 1.58 = 0.44303797468354 (passes)
past k_y = 1 the frequency falls: 0.4 at k_y = 2 (passes)

The growth rate is an eigenvalue

At finite $C$ the density lags the potential, and the lag lets the wave draw energy from the gradient. For one Fourier mode, $ϕ = - ζ / k^{2}$ and $\partial_{y} \to i k_{y}$ , and the linear equations become $\frac{d}{d t} (ζ, n) = M (ζ, n)$ with the $2 \times 2$ matrix under the title. Its eigenvalues are $λ = γ - iω$ : the real part is the growth rate, and minus the imaginary part is the frequency.

The trace of $M$ is real and its determinant, $i C κ k_{y} / k^{2}$ , is imaginary, so the quadratic formula needs the square root of a complex number. Viv has complex arithmetic, and sqrt of a complex value is the principal root, whose real part is not negative; that picks the eigenvalue that grows faster.

At $C = 1$ and $k = (0, 1)$ : $γ = 0.098684$ and $ω = 0.45509$ , below the drift-wave 0.5. Scanning $k_{y}$ , the fastest growth at $C = 1$ is $γ = 0.10687$ near $k_{y} = 1.30$ .

new: hw_lambda, growth, freq

Rung 2 tests 3 of 3 pass

at C = 1, k = (0, 1): gamma = 0.09868411347 and omega = 0.45508986056 (numpy eigvals) (passes)
at C = 0.5, k = (0.3, 0.7): gamma = 0.10155874337 and omega = 0.38554401937 (numpy eigvals) (passes)
the fastest growth at C = 1 is 0.106870 near k_y = 1.30 (numpy scan with step 0.0005) (passes)

Stiff coupling gives back the drift wave

Rung 1 took $C \to \infty$ by hand. The eigenvalue should agree, and expanding the small root of $λ^{2} - tr M λ + det M = 0$ in powers of $1/ C$ says how fast. With $tr M = - C (1 + k^{2}) / k^{2}$ , the first term is $det M / tr M = - i ω_{H M}$ , the drift wave. The next two terms are

$γ \approx \frac{ω _{H M}^{2} k ^{2}}{C ( 1 + k ^{2} )}, 1 - \frac{ω}{ω _{H M}} \approx \frac{2 ω _{H M}^{2} k ^{4}}{C ^{2} ( 1 + k ^{2} ) ^{2}} .$

The growth rate vanishes like $1/ C$ and the frequency error like $1/ C^{2}$ . At $k = (0, 1)$ both coefficients are $1/8$ , and the cell checks $C γ$ and $C^{2} (1 - ω / ω_{H M})$ at $C = 1000$ . The growth needs the lag: electrons that keep up exactly give a wave that neither grows nor decays.

Rung 3 tests 2 of 2 pass

C gamma approaches omega_HM^2 k^2 / (1 + k^2) = 1/8 (0.12499996 at C = 1000) (passes)
C^2 (1 - omega / omega_HM) approaches 2 omega_HM^2 k^4 / (1 + k^2)^2 = 1/8 (passes)

A zonal flow relaxes to k^2 / (1 + k^2) of its start

A zonal flow is a mode with $k_{y} = 0$ : a band of $E \times B$ flow in the $y$ direction, constant along $y$ . For it the drive $κ \partial_{y} ϕ$ is zero and only the coupling acts, $\partial_{t} ζ = \partial_{t} n = C (- ζ / k^{2} - n)$ . Here $M$ is singular. One eigenvalue is $- C (1 + 1/ k^{2})$ ; the other is 0, because $ζ - n$ never changes. So a band started with $ζ = 1$ , $n = 0$ does not decay away: it relaxes to

$\frac{ζ ( t )}{ζ ( 0 )} = \frac{k ^{2}}{1 + k ^{2}} + \frac{e ^{- C (1 + 1/ k^{2}) t}}{1 + k ^{2}},$

which is $0.2$ at $k_{x} = 0.5$ . Wide bands (small $k$ ) keep little, narrow ones keep most. The cell integrates the two equations with fourth-order Runge-Kutta to $t = 2$ , where the exact ratio is $0.2000363$ .

new: zonal_history

Rung 4 tests 2 of 2 pass

Runge-Kutta gives the exact k^2/(1+k^2) + e^(-C(1+1/k^2)t)/(1+k^2) at t = 2 to 1e-12 (passes)
and scipy's matrix exponential gives 0.20003631994381 for the same band (passes)

The modified model leaves the band alone

In a tokamak the coupling term comes from electrons moving along field lines, and a zonal flow is constant on a whole flux surface, so there is nothing along the field line for the electrons to smooth out. The modified Hasegawa-Wakatani model (Numata, Ball and Dewar, 2007) writes this in: the coupling acts only on the non-zonal parts, $C (ϕ ~ - n ~)$ , where the tilde removes the average over $y$ . For a $k_{y} = 0$ mode that makes the coupling exactly zero, and every non-zonal mode is unchanged.

With the same Runge-Kutta loop the zonal band now stays at its starting value, where rung 4's fell to $0.2$ . The band has no linear sink left; in the full equations it changes through the nonlinear bracket: the Reynolds stress of the drift waves around it.

new: modified_coupling

Rung 5 tests 2 of 2 pass

in the modified model the zonal band keeps its full amplitude (passes)
a non-zonal mode is unchanged: at k = (0, 1) the growth rate is still 0.09868411347 (passes)

Printed

w1: 0.5

0.443038

0.0986841

0.45509

g_best: 0.106869

ky_best: 1.3

c_gamma: 0.125

c2_dw: 0.125

zr: 0.200036

zr_exact: 0.200036

zr_mod: 1

classical Hasegawa-Wakatani

The same band in the two models

classicalmodified

1The drift wave

2The growth rate is an eigenvalue

3Stiff coupling gives back the drift wave

4A zonal flow relaxes to k^2 / (1 + k^2) of its start

5The modified model leaves the band alone

The drift wave

The growth rate is an eigenvalue

Stiff coupling gives back the drift wave

A zonal flow relaxes to k^2 / (1 + k^2) of its start

The modified model leaves the band alone