Landau Damping Simulator

1. Real‑world phenomenon

Landau damping describes how electrostatic waves in a collisionless plasma lose amplitude because resonant particles surf the wave and extract energy. Unlike viscous or resistive losses, no thermalization or particle collisions are required; the wave hands energy directly to particles traveling at its phase velocity. This subtle mechanism is pivotal in fusion devices, space plasmas, and accelerator beams where tiny electric fields shape macroscopic behavior. Observing Landau damping in the laboratory is difficult, so an interactive visualization can bridge the gap between formal theory and physical intuition. By animating a Langmuir wave whose amplitude decays according to the kinetic damping rate, the simulator lets users watch energy drain away in real time. The canvas makes the invisible electric oscillation tangible: users can adjust plasma parameters, press play, and observe how faster thermal speeds or higher wavenumbers accelerate the decay. Such interactivity turns an otherwise esoteric topic into a kinetic laboratory accessible from any browser.

2. Variables and assumptions

The model follows a homogeneous, unmagnetized, Maxwellian electron plasma with immobile ions. Input parameters include electron density $\mathrm{n\_e}$, electron temperature $\mathrm{T\_e}$, and wavenumber $k$. From these we compute the plasma frequency $\mathrm{\omega \_p}$, thermal velocity $\mathrm{v\_\{th\}}$, Debye length $\mathrm{\lambda \_D}$, and Landau damping rate $\gamma$. The wave is assumed electrostatic and longitudinal, so magnetic forces are absent. Amplitude is normalized to unity at t=0, and we track dimensionless energy proportional to the square of amplitude. Units are strictly SI: density in m⁻³, temperature in kelvin, wavenumber in m⁻¹, and time in seconds. Negative or non‑finite inputs are rejected. Because Landau damping is inherently linear, we restrict amplitudes to be small so that the wave does not modify the particle distribution. Spatial variation is represented by a single sinusoid—effects such as wave packet dispersion, nonlinear trapping, or multi‑dimensional structure are beyond the scope of this toy model.

3. Governing equations

The electron plasma frequency is $$\mathrm{\omega \_p}=\sqrt{\frac{\mathrm{n\_e}{e}^2}{\mathrm{\epsilon \_0}\mathrm{m\_e}}}$$. The thermal velocity is $\mathrm{v\_\{th\}}=\sqrt{\frac{\mathrm{k\_B}\mathrm{T\_e}}{\mathrm{m\_e}}}$, leading to the Debye length $\mathrm{\lambda \_D}=\sqrt{\frac{\mathrm{\epsilon \_0}\mathrm{k\_B}\mathrm{T\_e}}{\mathrm{n\_e}{e}^2}}$. The classic Landau damping rate for Langmuir waves with phase velocity exceeding the thermal speed is

$$\gamma =-\frac{\sqrt{\pi }}{2}\mathrm{\omega \_p}{k\mathrm{\lambda \_D}}^{-3}{e}^{-\frac{1}{2}{k\mathrm{\lambda \_D}}^2}$$

The wave frequency is approximated by $\omega =\sqrt{\mathrm{\omega \_p}{}^{}+\frac{3}{2}\frac{k^2}{{\mathrm{v\_\{th\}}}^2}}\mathrm{\omega \_p}$, sufficient for weak damping. The envelope amplitude $A$ obeys $\frac{\mathrm{dA}}{\mathrm{dt}}=\gamma A$, whose solution is exponential decay $\mathrm{A(t)}=\mathrm{A\_0}{e}^{\mathrm{\gamma t}}$. Energy density scales as $E\propto A^2$. The simulator integrates the amplitude equation numerically and displays both the spatial wave $\mathrm{\phi (x,t)=A(t)\backslash sin(kx-\omega t)}$ and the normalized energy $\mathrm{E/E\_0}$ via a striped bar.

4. Numerical scheme

An explicit Euler method advances the amplitude: $\mathrm{A\_\{n+1\}}=\mathrm{A\_n}+\gamma \mathrm{A\_n}\mathrm{\Delta t}$. The time step $\mathrm{\Delta t}$ must satisfy $\mathrm{|\gamma |\Delta t\; \backslash ll\; 1}$ for accuracy; the input is clamped between 10⁻⁷ and 10⁻³ s. The phase $\varphi$ of the sinusoid evolves via $\mathrm{\varphi \_\{n+1\}=\varphi \_n+\omega \Delta t}$. Although the amplitude equation admits an analytical solution, integrating numerically highlights stability issues: overly large steps exaggerate or even invert the decay. The script tracks initial energy and reports fractional drift between numerical and analytic amplitudes to alert users when $\mathrm{\Delta t}$ is too coarse. Because Euler’s global error scales with $\mathrm{\Delta t}$, halving the step roughly halves the energy discrepancy.

5. Worked example

Consider a laboratory plasma with $\mathrm{n\_e=10^\{18\}}$ m⁻³, $\mathrm{T\_e=10^5}$ K, and $\mathrm{k=1}$ m⁻¹. Plugging these values into the formulas yields $\mathrm{\omega \_p\approx 5.6\times 10^\{10\}}$ s⁻¹, $\mathrm{\lambda \_D\approx 2\times 10^\{-3\}}$ m, and a damping rate of $\mathrm{\gamma \approx -1.2\times 10^5}$ s⁻¹. With $\mathrm{\Delta t=10^\{-6\}}$ s, the numerical integration produces a wave that decays to half its initial amplitude after about 5.8 microseconds. The canvas plots the oscillation at each step, while the caption reports the current amplitude and theoretical half‑life. Exporting the CSV provides time‑series data of amplitude and energy for plotting in external software. The energy bar shrinks in sync with the exponential decay, reinforcing the link between amplitude and stored field energy.

6. Comparison table

The table compares the baseline parameters above with two variants: one warmer plasma and one shorter wavelength.

Scenario	Tₑ (K)	k (m⁻¹)	|γ| (s⁻¹)	Half‑life (μs)
Baseline	1×10⁵	1	1.2×10⁵	5.8
Hot plasma	5×10⁵	1	2.7×10⁵	2.6
Short wavelength	1×10⁵	5	4.8×10⁶	0.14

The warmer plasma increases the thermal velocity, steepening the distribution function at the phase velocity and thus boosting the damping rate. Raising the wavenumber has an even stronger effect because of the ${k\mathrm{\lambda \_D}}^{-3}$ dependence, causing the wave to vanish almost instantaneously.

7. How to read the animation

The canvas displays the wave as an orange curve oscillating about the midline. Time runs from left to right as the phase advances; amplitude shrinkage is visible as the peaks collapse toward the center. The striped bar beneath the canvas represents remaining wave energy relative to the start. Keyboard users can focus the canvas and press the space bar to toggle play and pause. The caption narrates current time, amplitude, and energy fraction, providing parity for screen‑reader users. When the amplitude falls below 1% of the initial value, the simulation stops and the caption summarizes the total energy lost.

8. Limitations

The simulator assumes a Maxwellian electron distribution and neglects magnetic fields, nonlinear wave‑particle trapping, relativistic effects, and spatial inhomogeneities. The exponential damping formula is valid only for small $\mathrm{k\lambda \_D}$ and weak damping. Large amplitudes or non‑Maxwellian tails can lead to growth or plateau formation that the model cannot capture. Numerically, explicit Euler requires small time steps; selecting a large $\mathrm{\Delta t}$ may produce artificial growth or negative energies. The visualized wave is one‑dimensional, so effects like oblique propagation or mode coupling are omitted.

9. Suggested extensions

Future versions could integrate the full Vlasov–Poisson system using a particle‑in‑cell scheme, animate phase‑space vortices, or include ion Landau damping for ion‑acoustic waves. Introducing external electric fields would allow exploration of inverse Landau damping and wave growth. A phase‑space plot overlay could display the resonant velocity region, while adaptive time stepping would improve efficiency when damping is slow. Because the code is entirely client‑side, enthusiasts can fork it to experiment with alternative distributions or magnetic geometries.

10. References and related tools

For detailed derivations, see Landau’s original paper (J. Phys. USSR 10, 25, 1946) or textbooks such as Nicholson’s Introduction to Plasma Theory. This simulator complements our Lorentz Force Calculator, the Mean Free Path Calculator, and the Point Charge Field Calculator, which explore other facets of plasma and electromagnetic physics.