1#!/usr/bin/env python3
2"""
3Landau Damping Simulation
4
5This script demonstrates Landau damping using the 1D Vlasov-Poisson solver.
6Shows exponential damping of electric field energy and phase space filamentation.
7Also demonstrates inverse Landau damping (bump-on-tail instability).
8
9Author: Plasma Physics Examples
10License: MIT
11"""
12
13import numpy as np
14import matplotlib.pyplot as plt
15from scipy.constants import e, m_e, epsilon_0, k
16from scipy.interpolate import CubicSpline
17from scipy.fft import fft, ifft, fftfreq
18from scipy.special import erf
19
20
21class VlasovPoisson1D:
22 """1D-1V Vlasov-Poisson solver (same as vlasov_1d.py)."""
23
24 def __init__(self, Nx, Nv, Lx, v_max, n0, T, q=-e, m=m_e):
25 self.Nx = Nx
26 self.Nv = Nv
27 self.Lx = Lx
28 self.v_max = v_max
29 self.n0 = n0
30 self.T = T
31 self.q = q
32 self.m = m
33
34 self.x = np.linspace(0, Lx, Nx, endpoint=False)
35 self.dx = Lx / Nx
36 self.v = np.linspace(-v_max, v_max, Nv)
37 self.dv = 2 * v_max / Nv
38
39 self.f = np.zeros((Nx, Nv))
40 self.E = np.zeros(Nx)
41 self.kx = 2 * np.pi * fftfreq(Nx, Lx / Nx)
42
43 def initialize_maxwellian(self, perturbation=0.0, k_pert=1):
44 v_th = np.sqrt(2 * k * self.T / self.m)
45 f_M = (self.n0 / (np.sqrt(2 * np.pi) * v_th) *
46 np.exp(-self.v**2 / (2 * v_th**2)))
47 k_wave = 2 * np.pi * k_pert / self.Lx
48 n_pert = 1 + perturbation * np.cos(k_wave * self.x)
49 for i in range(self.Nx):
50 self.f[i, :] = n_pert[i] * f_M
51
52 def initialize_bump_on_tail(self, n_beam_frac=0.1, v_beam_factor=3.0,
53 T_beam_factor=0.5, perturbation=0.01, k_pert=1):
54 """Initialize with bump-on-tail distribution."""
55 v_th = np.sqrt(2 * k * self.T / self.m)
56 v_beam = v_beam_factor * v_th
57 T_beam = T_beam_factor * self.T
58 v_th_beam = np.sqrt(2 * k * T_beam / self.m)
59
60 n_bg = self.n0 * (1 - n_beam_frac)
61 n_beam = self.n0 * n_beam_frac
62
63 # Background Maxwellian
64 f_bg = (n_bg / (np.sqrt(2 * np.pi) * v_th) *
65 np.exp(-self.v**2 / (2 * v_th**2)))
66
67 # Beam
68 f_beam = (n_beam / (np.sqrt(2 * np.pi) * v_th_beam) *
69 np.exp(-(self.v - v_beam)**2 / (2 * v_th_beam**2)))
70
71 # Total
72 f_total = f_bg + f_beam
73
74 # Spatial perturbation
75 k_wave = 2 * np.pi * k_pert / self.Lx
76 n_pert = 1 + perturbation * np.cos(k_wave * self.x)
77
78 for i in range(self.Nx):
79 self.f[i, :] = n_pert[i] * f_total
80
81 def compute_density(self):
82 return np.trapz(self.f, self.v, axis=1)
83
84 def compute_electric_field(self):
85 n = self.compute_density()
86 rho = self.q * (n - self.n0)
87 rho_k = fft(rho)
88 E_k = np.zeros_like(rho_k, dtype=complex)
89 E_k[1:] = -1j * rho_k[1:] / (self.kx[1:] * epsilon_0)
90 E_k[0] = 0
91 self.E = np.real(ifft(E_k))
92
93 def advect_x(self, dt):
94 f_new = np.zeros_like(self.f)
95 for j in range(self.Nv):
96 x_old = (self.x - self.v[j] * dt) % self.Lx
97 cs = CubicSpline(self.x, self.f[:, j], bc_type='periodic')
98 f_new[:, j] = cs(x_old)
99 self.f = f_new
100
101 def advect_v(self, dt):
102 f_new = np.zeros_like(self.f)
103 a = self.q * self.E / self.m
104 for i in range(self.Nx):
105 v_old = self.v - a[i] * dt
106 mask = (v_old >= -self.v_max) & (v_old <= self.v_max)
107 if np.any(mask):
108 cs = CubicSpline(self.v, self.f[i, :], bc_type='natural')
109 f_new[i, mask] = cs(v_old[mask])
110 f_new[i, mask] = np.maximum(f_new[i, mask], 0)
111 self.f = f_new
112
113 def step_strang(self, dt):
114 self.advect_x(dt / 2)
115 self.compute_electric_field()
116 self.advect_v(dt)
117 self.advect_x(dt / 2)
118
119 def run(self, t_end, dt, save_interval=10):
120 n_steps = int(t_end / dt)
121 t_history = np.zeros(n_steps + 1)
122 E_history = np.zeros((n_steps + 1, self.Nx))
123 f_history = []
124
125 self.compute_electric_field()
126 t_history[0] = 0
127 E_history[0, :] = self.E
128 f_history.append(self.f.copy())
129
130 for n in range(n_steps):
131 self.step_strang(dt)
132 t_history[n + 1] = (n + 1) * dt
133 E_history[n + 1, :] = self.E
134
135 if n % save_interval == 0:
136 f_history.append(self.f.copy())
137
138 return {
139 't': t_history,
140 'E': E_history,
141 'f_snapshots': f_history,
142 'x': self.x,
143 'v': self.v
144 }
145
146
147def landau_damping_rate_theory(k, v_th):
148 """
149 Theoretical Landau damping rate (small k limit).
150
151 γ_L ≈ -sqrt(π/8) * (ω_pe/k)³ * (1/v_th³) * exp(-1/(2k²λ_D²) - 3/2)
152
153 For k*λ_D << 1 (long wavelength).
154
155 Parameters:
156 -----------
157 k : float
158 Wavenumber [m^-1]
159 v_th : float
160 Thermal velocity [m/s]
161
162 Returns:
163 --------
164 gamma : float
165 Damping rate [s^-1]
166 """
167 # Using simplified formula for k*lambda_D ~ 0.5
168 # γ ≈ -sqrt(π/8) * ω_pe * exp(-1/(2k²λ_D²) - 3/2)
169
170 # For the standard case, use empirical fit
171 # This is approximate; exact value requires solving dispersion relation
172 omega_pe = np.sqrt(1e16 * e**2 / (epsilon_0 * m_e)) # Nominal
173 k_lambda_D = k * v_th / omega_pe # Approximation
174
175 gamma = -np.sqrt(np.pi / 8) * omega_pe * np.exp(-1 / (2 * k_lambda_D**2) - 1.5)
176
177 return gamma
178
179
180def simulate_landau_damping():
181 """Simulate Landau damping."""
182
183 # Parameters
184 n0 = 1e16 # m^-3
185 T = 1e4 # K
186 v_th = np.sqrt(2 * k * T / m_e)
187
188 omega_pe = np.sqrt(n0 * e**2 / (epsilon_0 * m_e))
189 T_pe = 2 * np.pi / omega_pe
190 lambda_D = np.sqrt(epsilon_0 * k * T / (n0 * e**2))
191
192 print(f"\nLandau Damping Simulation:")
193 print(f" Plasma frequency: {omega_pe/(2*np.pi):.2e} Hz")
194 print(f" Debye length: {lambda_D:.2e} m")
195 print(f" Thermal velocity: {v_th:.2e} m/s")
196
197 # Grid
198 Nx = 64
199 Nv = 256 # Higher resolution for damping
200 Lx = 10 * lambda_D
201 v_max = 6 * v_th
202
203 k_pert = 1
204 k_wave = 2 * np.pi * k_pert / Lx
205
206 # Theoretical damping rate
207 gamma_theory = landau_damping_rate_theory(k_wave, v_th)
208 print(f" Wavenumber k: {k_wave:.2e} m^-1")
209 print(f" k*λ_D: {k_wave * lambda_D:.3f}")
210 print(f" Theoretical damping rate: {gamma_theory:.2e} s^-1")
211 print(f" Damping time: {-1/gamma_theory:.2e} s ({-1/(gamma_theory*T_pe):.2f} T_pe)")
212
213 # Initialize solver
214 solver = VlasovPoisson1D(Nx, Nv, Lx, v_max, n0, T)
215 solver.initialize_maxwellian(perturbation=0.1, k_pert=k_pert)
216
217 # Run
218 t_end = 50 * T_pe
219 dt = T_pe / 50
220
221 print(f" Running for {t_end/T_pe:.1f} T_pe...")
222
223 history = solver.run(t_end, dt, save_interval=max(1, int(t_end/dt) // 50))
224
225 return history, omega_pe, lambda_D, gamma_theory
226
227
228def simulate_bump_on_tail():
229 """Simulate bump-on-tail instability (inverse Landau damping)."""
230
231 # Parameters
232 n0 = 1e16
233 T = 1e4
234 v_th = np.sqrt(2 * k * T / m_e)
235
236 omega_pe = np.sqrt(n0 * e**2 / (epsilon_0 * m_e))
237 T_pe = 2 * np.pi / omega_pe
238 lambda_D = np.sqrt(epsilon_0 * k * T / (n0 * e**2))
239
240 print(f"\nBump-on-Tail Simulation:")
241 print(f" Beam fraction: 10%")
242 print(f" Beam velocity: 3 v_th")
243
244 # Grid
245 Nx = 64
246 Nv = 256
247 Lx = 10 * lambda_D
248 v_max = 8 * v_th # Wider for beam
249
250 # Initialize with bump-on-tail
251 solver = VlasovPoisson1D(Nx, Nv, Lx, v_max, n0, T)
252 solver.initialize_bump_on_tail(n_beam_frac=0.1, v_beam_factor=3.0,
253 T_beam_factor=0.5, perturbation=0.01, k_pert=1)
254
255 # Run
256 t_end = 100 * T_pe
257 dt = T_pe / 50
258
259 print(f" Running for {t_end/T_pe:.1f} T_pe...")
260
261 history = solver.run(t_end, dt, save_interval=max(1, int(t_end/dt) // 50))
262
263 return history, omega_pe, lambda_D
264
265
266def plot_landau_damping(history, omega_pe, lambda_D, gamma_theory):
267 """Plot Landau damping results."""
268
269 t = history['t']
270 E = history['E']
271 x = history['x']
272 v = history['v']
273 f_snapshots = history['f_snapshots']
274
275 T_pe = 2 * np.pi / omega_pe
276 v_th = np.sqrt(2 * k * 1e4 / m_e)
277
278 # Electric field energy
279 E_energy = np.sum(E**2, axis=1)
280
281 # Fit exponential decay
282 log_E = np.log(E_energy + 1e-20)
283 t_fit = t[t < 30 * T_pe]
284 log_E_fit = log_E[t < 30 * T_pe]
285
286 # Linear fit to log(E)
287 p = np.polyfit(t_fit, log_E_fit, 1)
288 gamma_measured = p[0]
289
290 print(f"\nMeasured damping rate: {gamma_measured:.2e} s^-1")
291 print(f"Theory: {gamma_theory:.2e} s^-1")
292 print(f"Error: {abs(gamma_measured - gamma_theory) / abs(gamma_theory) * 100:.1f}%")
293
294 # Create figure
295 fig = plt.figure(figsize=(16, 10))
296
297 # Plot 1: E field energy (log scale)
298 ax1 = plt.subplot(2, 3, 1)
299 ax1.semilogy(t / T_pe, E_energy, 'b-', linewidth=2, label='Simulation')
300 ax1.semilogy(t_fit / T_pe, np.exp(p[0] * t_fit + p[1]), 'r--', linewidth=2,
301 label=f'Fit: γ = {gamma_measured:.2e} s⁻¹')
302 ax1.semilogy(t / T_pe, E_energy[0] * np.exp(2 * gamma_theory * t), 'g:',
303 linewidth=2, label=f'Theory: γ = {gamma_theory:.2e} s⁻¹')
304
305 ax1.set_xlabel('Time [T_pe]', fontsize=11)
306 ax1.set_ylabel('∫E² dx [a.u.]', fontsize=11)
307 ax1.set_title('Electric Field Energy Decay', fontsize=12, fontweight='bold')
308 ax1.legend(loc='best', fontsize=9)
309 ax1.grid(True, alpha=0.3)
310
311 # Plot 2-7: Phase space evolution
312 snapshot_times = [0, 10, 20, 30, 40, 50]
313
314 for idx, t_snap_idx in enumerate(range(min(6, len(f_snapshots)))):
315 if t_snap_idx >= len(f_snapshots):
316 break
317
318 ax = plt.subplot(2, 3, 2 + idx)
319 f = f_snapshots[t_snap_idx]
320 t_actual = t_snap_idx * (t[-1] / (len(f_snapshots) - 1))
321
322 cs = ax.contourf(x / lambda_D, v / v_th, f.T, levels=30, cmap='hot')
323 ax.set_xlabel('x [λ_D]', fontsize=10)
324 ax.set_ylabel('v [v_th]', fontsize=10)
325 ax.set_title(f't = {t_actual/T_pe:.1f} T_pe', fontsize=11, fontweight='bold')
326 ax.set_ylim([-4, 4])
327
328 plt.tight_layout()
329 plt.savefig('landau_damping.png', dpi=300, bbox_inches='tight')
330 plt.show()
331
332
333def plot_bump_on_tail(history, omega_pe, lambda_D):
334 """Plot bump-on-tail (growth) results."""
335
336 t = history['t']
337 E = history['E']
338 x = history['x']
339 v = history['v']
340 f_snapshots = history['f_snapshots']
341
342 T_pe = 2 * np.pi / omega_pe
343 v_th = np.sqrt(2 * k * 1e4 / m_e)
344
345 # Electric field energy
346 E_energy = np.sum(E**2, axis=1)
347
348 # Fit exponential growth (linear phase)
349 t_linear = (t > 10 * T_pe) & (t < 40 * T_pe)
350 log_E = np.log(E_energy[t_linear] + 1e-20)
351 p = np.polyfit(t[t_linear], log_E, 1)
352 gamma_growth = p[0]
353
354 print(f"\nMeasured growth rate: {gamma_growth:.2e} s^-1")
355
356 # Create figure
357 fig = plt.figure(figsize=(16, 10))
358
359 # Plot 1: E field energy (log scale, showing growth)
360 ax1 = plt.subplot(2, 3, 1)
361 ax1.semilogy(t / T_pe, E_energy, 'b-', linewidth=2, label='Simulation')
362 ax1.semilogy(t[t_linear] / T_pe, np.exp(p[0] * t[t_linear] + p[1]), 'r--',
363 linewidth=2, label=f'Growth: γ = {gamma_growth:.2e} s⁻¹')
364
365 ax1.set_xlabel('Time [T_pe]', fontsize=11)
366 ax1.set_ylabel('∫E² dx [a.u.]', fontsize=11)
367 ax1.set_title('Electric Field Energy Growth (Instability)', fontsize=12, fontweight='bold')
368 ax1.legend(loc='best', fontsize=10)
369 ax1.grid(True, alpha=0.3)
370
371 # Plot 2-7: Phase space showing plateau formation
372 for idx in range(min(6, len(f_snapshots))):
373 ax = plt.subplot(2, 3, 2 + idx)
374 f = f_snapshots[idx]
375 t_actual = idx * (t[-1] / (len(f_snapshots) - 1))
376
377 cs = ax.contourf(x / lambda_D, v / v_th, f.T, levels=30, cmap='hot')
378 ax.set_xlabel('x [λ_D]', fontsize=10)
379 ax.set_ylabel('v [v_th]', fontsize=10)
380 ax.set_title(f't = {t_actual/T_pe:.1f} T_pe', fontsize=11, fontweight='bold')
381 ax.set_ylim([-5, 6])
382
383 plt.tight_layout()
384 plt.savefig('bump_on_tail_growth.png', dpi=300, bbox_inches='tight')
385 plt.show()
386
387
388if __name__ == '__main__':
389 print("\n" + "="*80)
390 print("LANDAU DAMPING AND BUMP-ON-TAIL INSTABILITY")
391 print("="*80)
392
393 print("\nPart 1: Landau Damping (Collisionless Damping)")
394 print("-" * 80)
395 history_ld, omega_pe, lambda_D, gamma_theory = simulate_landau_damping()
396 plot_landau_damping(history_ld, omega_pe, lambda_D, gamma_theory)
397
398 print("\nPart 2: Bump-on-Tail (Inverse Landau Damping)")
399 print("-" * 80)
400 history_bt, omega_pe, lambda_D = simulate_bump_on_tail()
401 plot_bump_on_tail(history_bt, omega_pe, lambda_D)
402
403 print("\nDone! Generated 2 figures:")
404 print(" - landau_damping.png")
405 print(" - bump_on_tail_growth.png")