Part 8, Chapter 2

Gyrokinetics

Reduced kinetic theory for magnetized plasmas

2.1 Introduction and Motivation

The full six-dimensional kinetic equation is often computationally intractable for strongly magnetized plasmas. In tokamaks, particles gyrate rapidly around field lines at the cyclotron frequency, but the turbulence that drives transport evolves on much slower timescales. Gyrokinetics exploits this scale separation by averaging out the fast gyromotion, reducing the problem from six to five phase-space dimensions.

The key physical insight is that in a strongly magnetized plasma, a charged particle's trajectory can be decomposed into a fast circular gyration around a guiding center and a slow drift of that guiding center. If perturbations vary slowly compared to the gyrofrequency, we can average over the gyrophase and still capture all the essential physics of drift-wave turbulence, microinstabilities, and anomalous transport.

2.2 Gyrokinetic Ordering

Gyrokinetics is built on a formal ordering in a small parameter. We define:

$$\epsilon_{GK} \sim \frac{\omega}{\Omega_i} \sim \frac{k_\parallel}{k_\perp} \sim \frac{\rho_i}{L} \sim \frac{e\phi}{T_i} \ll 1$$

where $\Omega_i = eB/m_i$ is the ion cyclotron frequency,$\rho_i = v_{ti}/\Omega_i$ is the ion Larmor radius, and$L$ is the equilibrium gradient scale length.

The crucial feature is that perpendicular wavelengths are comparable to the Larmor radius:

$$k_\perp \rho_i \sim 1$$

while the frequency is much smaller than the cyclotron frequency:

$$\frac{\omega}{\Omega_i} \ll 1$$

This ordering captures drift-wave turbulence, which has fine perpendicular structure (comparable to the Larmor radius) but elongated parallel structure along field lines. The fluctuation amplitudes are small compared to equilibrium quantities, so we work perturbatively, but finite Larmor radius (FLR) effects must be retained to all orders in$k_\perp \rho_i$.

2.3 The Gyroaveraging Operator

The central mathematical tool in gyrokinetics is the gyroaveraging operator. For any field quantity evaluated at the particle position $\mathbf{x} = \mathbf{R} + \boldsymbol{\rho}$, where $\mathbf{R}$ is the guiding center and$\boldsymbol{\rho}$ is the Larmor radius vector, the gyroaverage is:

$$\langle \phi \rangle_\alpha = \frac{1}{2\pi} \oint \phi(\mathbf{R} + \boldsymbol{\rho})\, d\alpha$$

where $\alpha$ is the gyrophase angle. In Fourier space, the gyroaveraging becomes multiplication by a Bessel function:

$$\langle \phi_\mathbf{k} \rangle_\alpha = J_0(k_\perp \rho_i)\, \phi_\mathbf{k}$$

The Bessel function $J_0(k_\perp \rho_i)$ encodes all finite Larmor radius effects. For long wavelengths ($k_\perp\rho_i \ll 1$),$J_0 \approx 1 - k_\perp^2\rho_i^2/4$, recovering the drift-kinetic limit. For short wavelengths ($k_\perp\rho_i \gg 1$), the gyroaverage strongly suppresses the fluctuations as the particle samples many oscillations around its orbit.

2.4 The Gyrokinetic Vlasov Equation

After performing the gyrophase average and a near-identity (Lie) coordinate transformation to remove fast oscillations order by order, the distribution function$f(\mathbf{R}, v_\parallel, \mu, t)$ evolves according to the gyrokinetic equation:

$$\frac{\partial f}{\partial t} + \dot{\mathbf{R}} \cdot \nabla f + \dot{v}_\parallel \frac{\partial f}{\partial v_\parallel} = C[f]$$

The guiding center velocity contains the parallel streaming, the $E \times B$ drift, and curvature/grad-B drifts:

$$\dot{\mathbf{R}} = v_\parallel \hat{\mathbf{b}} + \frac{\hat{\mathbf{b}}}{B\Omega_s} \times \left(\mu \nabla B + v_\parallel^2 \hat{\mathbf{b}} \cdot \nabla \hat{\mathbf{b}} + \frac{q_s}{m_s} \nabla \langle \phi \rangle_\alpha \right)$$

The parallel acceleration is:

$$\dot{v}_\parallel = -\frac{1}{m_s}\left(\mu \nabla_\parallel B + q_s \nabla_\parallel \langle \phi \rangle_\alpha \right)$$

The magnetic moment $\mu = m v_\perp^2 / (2B)$ is an adiabatic invariant and does not evolve (to lowest order). This reduces the phase space from 6D to 5D: three guiding-center coordinates $\mathbf{R}$, the parallel velocity$v_\parallel$, and the magnetic moment$\mu$ (which enters as a parameter).

2.5 Ion Temperature Gradient (ITG) Mode

The ion temperature gradient (ITG) instability is a primary driver of anomalous transport in tokamaks. We derive it from the gyrokinetic equation by linearizing around a Maxwellian equilibrium with density and temperature gradients.

Consider a slab geometry with gradients in the $x$-direction and a uniform magnetic field $\mathbf{B} = B\hat{\mathbf{z}}$. The equilibrium distribution is:

$$F_0 = \frac{n_0(x)}{(2\pi T_i(x)/m_i)^{3/2}} \exp\!\left(-\frac{m_i v^2}{2T_i(x)}\right)$$

We define the gradient scale lengths:

$$\frac{1}{L_n} = -\frac{1}{n_0}\frac{dn_0}{dx}, \qquad \frac{1}{L_T} = -\frac{1}{T_i}\frac{dT_i}{dx}, \qquad \eta_i = \frac{L_n}{L_T}$$

The diamagnetic drift frequency and its gradient correction are:

$$\omega_{*i} = \frac{k_y T_i}{eB L_n}, \qquad \omega_{*Ti} = \omega_{*i}\,\eta_i$$

After linearizing the gyrokinetic equation and taking appropriate limits (flat density for the pure ITG, adiabatic electrons), the local dispersion relation becomes:

$$1 + \frac{1}{\tau} - \frac{\omega_{*i}}{\omega}\left[1 + \eta_i\left(\frac{\omega}{\omega_{*i}} \frac{v_{ti}^2 k_\parallel^2}{\omega^2} - \frac{1}{2}\right)\right]\Gamma_0(b) = 0$$

where $\tau = T_i/T_e$,$b = k_\perp^2 \rho_i^2$, and$\Gamma_0(b) = I_0(b)e^{-b}$ with $I_0$ the modified Bessel function. The ITG mode becomes unstable when $\eta_i$ exceeds a critical threshold, typically $\eta_{i,\text{crit}} \approx 2/3$ in slab geometry.

2.6 Gyrokinetic Field Equations

The gyrokinetic system is closed by the quasineutrality (gyrokinetic Poisson) equation. Since ions are gyroaveraged but electrons respond adiabatically at long wavelengths:

$$\frac{e n_0}{T_e}\left(\phi - \langle\phi\rangle_\text{fs}\right) = \int J_0(k_\perp \rho_i)\, \delta f_i\, d^3v$$

The double gyroaveraging on the right-hand side (once from the charge density, once from the field felt by the particles) produces the polarization density, leading to:

$$\frac{e^2 n_0}{T_i}\left(1 - \Gamma_0(b)\right)\phi = e \int J_0\, \delta h_i\, d^3v - \frac{e n_0}{T_e}\left(\phi - \langle\phi\rangle_\text{fs}\right)$$

Here $\delta h_i$ is the non-adiabatic part of the ion distribution. The term$1 - \Gamma_0(b)$ represents the polarization shielding from the difference between the guiding-center density and the particle density. This is analogous to the classical polarization drift but retains all orders in FLR effects.

Interactive Simulation

ITG Growth Rate vs eta_i

Python

script.py90 lines

import numpy as np
from scipy.special import i0e  # I_0(b)*exp(-b) = Gamma_0(b)

# ITG dispersion relation (simplified slab model)
# We solve for complex omega given eta_i
# Dispersion: omega^2 - omega * omega_star_i * (1 + eta_i*(v^2/v_ti^2 - 3/2)) * Gamma0 = k_par^2 * v_ti^2 * tau

# Parameters
T_i_eV = 1000.0      # Ion temperature in eV
T_e_eV = 1000.0      # Electron temperature in eV
tau = T_i_eV / T_e_eV
B0 = 3.0              # Tesla
m_i = 1.67e-27 * 2    # Deuterium mass
e = 1.602e-19
v_ti = np.sqrt(T_i_eV * e / m_i)
Omega_i = e * B0 / m_i
rho_i = v_ti / Omega_i

# Mode parameters
k_perp_rho = 0.5      # k_perp * rho_i
b = k_perp_rho**2
Gamma0 = i0e(b)       # I_0(b)*exp(-b)
L_n = 0.5             # Density gradient scale length (m)
k_y = k_perp_rho / rho_i
omega_star = k_y * T_i_eV * e / (e * B0 * L_n)  # Diamagnetic freq
k_par = 0.1 / 1.0     # Parallel wavevector (1/m)

# Scan eta_i
eta_i_vals = np.linspace(0.0, 6.0, 80)
gamma_vals = []
omega_r_vals = []

for eta_i in eta_i_vals:
    # Quadratic dispersion: omega^2 - omega * omega_star*(1 + eta_i*(E_avg - 3/2))*Gamma0
    #                        - k_par^2 * v_ti^2 * (1 + tau*eta_i) = 0
    # Simplified fluid ITG:
    # omega^2 + omega * omega_star * Gamma0 - omega * omega_star * eta_i * Gamma0
    #   - k_par^2 * v_ti^2 = 0
    # i.e. omega^2 - omega * omega_star * (eta_i - 1) * Gamma0 - k_par^2 * v_ti^2 = 0

A = 1.0
    B_coeff = -omega_star * (eta_i - 1.0) * Gamma0
    C = -k_par**2 * v_ti**2

disc = B_coeff**2 - 4*A*C
    if disc < 0:
        omega_plus = (-B_coeff + 1j*np.sqrt(-disc)) / (2*A)
    else:
        omega_plus = (-B_coeff + np.sqrt(disc)) / (2*A)

# For ITG: the unstable branch
    omega_complex = (-B_coeff - np.sqrt(disc + 0j)) / (2*A)
    gamma_vals.append(np.imag(omega_complex) / omega_star)
    omega_r_vals.append(np.real(omega_complex) / omega_star)

gamma_arr = np.array(gamma_vals)
omega_r_arr = np.array(omega_r_vals)

# Find threshold
unstable_mask = gamma_arr > 0.01
if np.any(unstable_mask):
    eta_crit = eta_i_vals[unstable_mask][0]
else:
    eta_crit = float('nan')

print(f"=== ITG Growth Rate Analysis ===")
print(f"k_perp * rho_i = {k_perp_rho:.2f}")
print(f"Gamma_0(b)     = {Gamma0:.4f}")
print(f"omega_*i       = {omega_star:.2e} rad/s")
print(f"k_parallel     = {k_par:.2f} /m")
print(f"v_ti           = {v_ti:.2e} m/s")
print(f"rho_i          = {rho_i*100:.2f} cm")
print(f"")
print(f"Critical eta_i ~ {eta_crit:.2f}")
print(f"")
print(f"{'eta_i':>8} {'gamma/omega*':>14} {'omega_r/omega*':>16}")
print("-" * 42)
for i in range(0, len(eta_i_vals), 8):
    print(f"{eta_i_vals[i]:8.2f} {gamma_arr[i]:14.4f} {omega_r_arr[i]:16.4f}")
print(f"")
print(f"Max growth rate: gamma/omega* = {np.max(gamma_arr):.4f} at eta_i = {eta_i_vals[np.argmax(gamma_arr)]:.2f}")

# Chart data
print(f"")
print(f"CHART_DATA:ITG Growth Rate vs eta_i")
print(f"x_label:eta_i = L_n / L_T")
print(f"y_label:gamma / omega_*i")
for ei, gi in zip(eta_i_vals, gamma_arr):
    print(f"POINT:{ei:.3f},{gi:.6f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Previous Next →