1#!/usr/bin/env python3
2"""
3Bump-on-Tail Instability - Detailed Analysis
4
5This script provides a detailed analysis of bump-on-tail instability including:
6- Linear growth phase
7- Nonlinear saturation
8- Plateau formation (quasilinear theory)
9- Particle trapping in wave potential
10
11Author: Plasma Physics Examples
12License: MIT
13"""
14
15import numpy as np
16import matplotlib.pyplot as plt
17from scipy.constants import e, m_e, epsilon_0, k
18from scipy.interpolate import CubicSpline
19from scipy.fft import fft, ifft, fftfreq
20
21
22class VlasovPoisson1D:
23 """1D-1V Vlasov-Poisson solver."""
24
25 def __init__(self, Nx, Nv, Lx, v_max, n0, T, q=-e, m=m_e):
26 self.Nx = Nx
27 self.Nv = Nv
28 self.Lx = Lx
29 self.v_max = v_max
30 self.n0 = n0
31 self.T = T
32 self.q = q
33 self.m = m
34
35 self.x = np.linspace(0, Lx, Nx, endpoint=False)
36 self.dx = Lx / Nx
37 self.v = np.linspace(-v_max, v_max, Nv)
38 self.dv = 2 * v_max / Nv
39
40 self.f = np.zeros((Nx, Nv))
41 self.E = np.zeros(Nx)
42 self.kx = 2 * np.pi * fftfreq(Nx, Lx / Nx)
43
44 def initialize_bump_on_tail(self, n_beam_frac=0.1, v_beam_factor=3.0,
45 T_beam_factor=0.5, perturbation=0.01, k_pert=1):
46 """Initialize with bump-on-tail distribution."""
47 v_th = np.sqrt(2 * k * self.T / self.m)
48 v_beam = v_beam_factor * v_th
49 T_beam = T_beam_factor * self.T
50 v_th_beam = np.sqrt(2 * k * T_beam / self.m)
51
52 n_bg = self.n0 * (1 - n_beam_frac)
53 n_beam = self.n0 * n_beam_frac
54
55 # Background
56 f_bg = (n_bg / (np.sqrt(2 * np.pi) * v_th) *
57 np.exp(-self.v**2 / (2 * v_th**2)))
58
59 # Beam
60 f_beam = (n_beam / (np.sqrt(2 * np.pi) * v_th_beam) *
61 np.exp(-(self.v - v_beam)**2 / (2 * v_th_beam**2)))
62
63 # Total
64 f_total = f_bg + f_beam
65
66 # Spatial perturbation
67 k_wave = 2 * np.pi * k_pert / self.Lx
68 n_pert = 1 + perturbation * np.cos(k_wave * self.x)
69
70 for i in range(self.Nx):
71 self.f[i, :] = n_pert[i] * f_total
72
73 # Store initial distribution for comparison
74 self.f_initial_v = f_total.copy()
75
76 def compute_density(self):
77 return np.trapz(self.f, self.v, axis=1)
78
79 def compute_electric_field(self):
80 n = self.compute_density()
81 rho = self.q * (n - self.n0)
82 rho_k = fft(rho)
83 E_k = np.zeros_like(rho_k, dtype=complex)
84 E_k[1:] = -1j * rho_k[1:] / (self.kx[1:] * epsilon_0)
85 E_k[0] = 0
86 self.E = np.real(ifft(E_k))
87
88 def compute_potential(self):
89 """Compute electrostatic potential φ from E = -dφ/dx."""
90 E_k = fft(self.E)
91 phi_k = np.zeros_like(E_k, dtype=complex)
92 phi_k[1:] = 1j * E_k[1:] / self.kx[1:]
93 phi_k[0] = 0
94 return np.real(ifft(phi_k))
95
96 def advect_x(self, dt):
97 f_new = np.zeros_like(self.f)
98 for j in range(self.Nv):
99 x_old = (self.x - self.v[j] * dt) % self.Lx
100 cs = CubicSpline(self.x, self.f[:, j], bc_type='periodic')
101 f_new[:, j] = cs(x_old)
102 self.f = f_new
103
104 def advect_v(self, dt):
105 f_new = np.zeros_like(self.f)
106 a = self.q * self.E / self.m
107 for i in range(self.Nx):
108 v_old = self.v - a[i] * dt
109 mask = (v_old >= -self.v_max) & (v_old <= self.v_max)
110 if np.any(mask):
111 cs = CubicSpline(self.v, self.f[i, :], bc_type='natural')
112 f_new[i, mask] = cs(v_old[mask])
113 f_new[i, mask] = np.maximum(f_new[i, mask], 0)
114 self.f = f_new
115
116 def step_strang(self, dt):
117 self.advect_x(dt / 2)
118 self.compute_electric_field()
119 self.advect_v(dt)
120 self.advect_x(dt / 2)
121
122 def run(self, t_end, dt, save_interval=10):
123 n_steps = int(t_end / dt)
124 t_history = np.zeros(n_steps + 1)
125 E_history = np.zeros((n_steps + 1, self.Nx))
126 phi_history = np.zeros((n_steps + 1, self.Nx))
127 f_history = []
128 f_avg_history = [] # Velocity-averaged f(v)
129
130 self.compute_electric_field()
131 t_history[0] = 0
132 E_history[0, :] = self.E
133 phi_history[0, :] = self.compute_potential()
134 f_history.append(self.f.copy())
135 f_avg_history.append(np.mean(self.f, axis=0))
136
137 for n in range(n_steps):
138 self.step_strang(dt)
139 t_history[n + 1] = (n + 1) * dt
140 E_history[n + 1, :] = self.E
141 phi_history[n + 1, :] = self.compute_potential()
142
143 if n % save_interval == 0:
144 f_history.append(self.f.copy())
145 f_avg_history.append(np.mean(self.f, axis=0))
146
147 return {
148 't': t_history,
149 'E': E_history,
150 'phi': phi_history,
151 'f_snapshots': f_history,
152 'f_avg': np.array(f_avg_history),
153 'x': self.x,
154 'v': self.v
155 }
156
157
158def simulate_full_evolution():
159 """Simulate full nonlinear evolution of bump-on-tail."""
160
161 # Parameters
162 n0 = 1e16
163 T = 1e4
164 v_th = np.sqrt(2 * k * T / m_e)
165
166 omega_pe = np.sqrt(n0 * e**2 / (epsilon_0 * m_e))
167 T_pe = 2 * np.pi / omega_pe
168 lambda_D = np.sqrt(epsilon_0 * k * T / (n0 * e**2))
169
170 print(f"\nBump-on-Tail Full Evolution:")
171 print(f" Thermal velocity: {v_th:.2e} m/s")
172 print(f" Plasma period: {T_pe:.2e} s")
173
174 # Grid
175 Nx = 128
176 Nv = 512 # High resolution for plateau formation
177 Lx = 20 * lambda_D
178 v_max = 8 * v_th
179
180 # Initialize
181 solver = VlasovPoisson1D(Nx, Nv, Lx, v_max, n0, T)
182 solver.initialize_bump_on_tail(n_beam_frac=0.15, v_beam_factor=3.0,
183 T_beam_factor=0.3, perturbation=0.02, k_pert=1)
184
185 # Run long enough to see saturation
186 t_end = 200 * T_pe
187 dt = T_pe / 50
188
189 print(f" Running for {t_end/T_pe:.1f} T_pe...")
190
191 history = solver.run(t_end, dt, save_interval=max(1, int(t_end/dt) // 100))
192
193 return history, solver, omega_pe, lambda_D, v_th
194
195
196def plot_full_evolution(history, solver, omega_pe, lambda_D, v_th):
197 """Plot complete evolution: linear, nonlinear, saturation."""
198
199 t = history['t']
200 E = history['E']
201 phi = history['phi']
202 x = history['x']
203 v = history['v']
204 f_snapshots = history['f_snapshots']
205 f_avg = history['f_avg']
206
207 T_pe = 2 * np.pi / omega_pe
208
209 # Electric field energy
210 E_energy = np.sum(E**2, axis=1)
211
212 # Create comprehensive figure
213 fig = plt.figure(figsize=(18, 14))
214
215 # Plot 1: Energy evolution (3 phases)
216 ax1 = plt.subplot(3, 4, 1)
217 ax1.semilogy(t / T_pe, E_energy, 'b-', linewidth=2)
218
219 # Mark phases
220 t_linear_end = 50
221 t_nonlinear_end = 120
222
223 ax1.axvspan(0, t_linear_end, alpha=0.2, color='green', label='Linear growth')
224 ax1.axvspan(t_linear_end, t_nonlinear_end, alpha=0.2, color='orange',
225 label='Nonlinear')
226 ax1.axvspan(t_nonlinear_end, t[-1]/T_pe, alpha=0.2, color='red',
227 label='Saturation')
228
229 ax1.set_xlabel('Time [T_pe]', fontsize=11)
230 ax1.set_ylabel('∫E² dx [a.u.]', fontsize=11)
231 ax1.set_title('Energy Evolution', fontsize=12, fontweight='bold')
232 ax1.legend(loc='best', fontsize=9)
233 ax1.grid(True, alpha=0.3)
234
235 # Plot 2: Distribution evolution in velocity
236 ax2 = plt.subplot(3, 4, 2)
237
238 # Plot initial, intermediate, and final distributions
239 time_indices = [0, len(f_avg)//3, 2*len(f_avg)//3, -1]
240 colors = ['blue', 'green', 'orange', 'red']
241 labels = ['Initial', 'Linear', 'Nonlinear', 'Saturated']
242
243 for idx, color, label in zip(time_indices, colors, labels):
244 ax2.plot(v / v_th, f_avg[idx] * v_th, color=color, linewidth=2,
245 label=label, alpha=0.8)
246
247 ax2.set_xlabel('v [v_th]', fontsize=11)
248 ax2.set_ylabel('f(v) × v_th', fontsize=11)
249 ax2.set_title('Distribution Function Evolution', fontsize=12, fontweight='bold')
250 ax2.legend(loc='best', fontsize=9)
251 ax2.grid(True, alpha=0.3)
252 ax2.set_xlim([-2, 6])
253
254 # Plot 3: df/dv showing slope flattening
255 ax3 = plt.subplot(3, 4, 3)
256
257 for idx, color, label in zip(time_indices, colors, labels):
258 dfdv = np.gradient(f_avg[idx], v)
259 ax3.plot(v / v_th, dfdv * v_th**2, color=color, linewidth=2,
260 label=label, alpha=0.8)
261
262 ax3.axhline(0, color='black', linestyle='--', linewidth=1)
263 ax3.set_xlabel('v [v_th]', fontsize=11)
264 ax3.set_ylabel('df/dv × v_th²', fontsize=11)
265 ax3.set_title('Slope (Plateau Formation)', fontsize=12, fontweight='bold')
266 ax3.legend(loc='best', fontsize=8)
267 ax3.grid(True, alpha=0.3)
268 ax3.set_xlim([1, 5])
269
270 # Plot 4: Growth rate measurement
271 ax4 = plt.subplot(3, 4, 4)
272
273 # Calculate instantaneous growth rate
274 log_E = np.log(E_energy + 1e-30)
275 gamma_inst = np.gradient(log_E, t)
276
277 ax4.plot(t / T_pe, gamma_inst * T_pe, 'b-', linewidth=2)
278 ax4.axhline(0, color='red', linestyle='--', linewidth=2)
279 ax4.set_xlabel('Time [T_pe]', fontsize=11)
280 ax4.set_ylabel('Growth Rate γ × T_pe', fontsize=11)
281 ax4.set_title('Instantaneous Growth Rate', fontsize=12, fontweight='bold')
282 ax4.grid(True, alpha=0.3)
283 ax4.set_ylim([-0.1, 0.3])
284
285 # Plots 5-12: Phase space snapshots at key times
286 snapshot_times = [0, 30, 60, 90, 120, 150, 180, -1]
287
288 for plot_idx, snap_idx in enumerate(snapshot_times):
289 ax = plt.subplot(3, 4, 5 + plot_idx)
290
291 if snap_idx == -1:
292 snap_idx = len(f_snapshots) - 1
293
294 if snap_idx < len(f_snapshots):
295 f = f_snapshots[snap_idx]
296 t_actual = snap_idx * (t[-1] / (len(f_snapshots) - 1))
297
298 # Determine phase
299 if t_actual < t_linear_end * T_pe:
300 phase = 'Linear'
301 elif t_actual < t_nonlinear_end * T_pe:
302 phase = 'Nonlinear'
303 else:
304 phase = 'Saturated'
305
306 cs = ax.contourf(x / lambda_D, v / v_th, f.T, levels=30, cmap='hot')
307 ax.set_xlabel('x [λ_D]', fontsize=9)
308 ax.set_ylabel('v [v_th]', fontsize=9)
309 ax.set_title(f'{phase}: t={t_actual/T_pe:.1f} T_pe',
310 fontsize=10, fontweight='bold')
311 ax.set_ylim([-3, 6])
312
313 plt.tight_layout()
314 plt.savefig('bump_on_tail_full_evolution.png', dpi=300, bbox_inches='tight')
315 plt.show()
316
317
318def plot_trapping_analysis(history, omega_pe, lambda_D, v_th):
319 """Analyze particle trapping in wave potential."""
320
321 t = history['t']
322 E = history['E']
323 phi = history['phi']
324 x = history['x']
325 v = history['v']
326 f_snapshots = history['f_snapshots']
327
328 T_pe = 2 * np.pi / omega_pe
329
330 # Select time in nonlinear phase (when trapping is strong)
331 t_trap_idx = len(f_snapshots) // 2
332
333 f_trap = f_snapshots[t_trap_idx]
334 phi_trap = phi[t_trap_idx * (len(t) // len(f_snapshots))]
335 E_trap = E[t_trap_idx * (len(t) // len(f_snapshots))]
336
337 # Wave phase velocity (estimate from k and omega)
338 k_wave = 2 * np.pi / (x[-1] - x[0])
339 # Approximate omega from df/dv=0 location (resonant velocity)
340 v_phase = 3 * v_th # Approximate
341
342 # Trapping width
343 phi_amplitude = (np.max(phi_trap) - np.min(phi_trap)) / 2
344 omega_bounce = np.sqrt(abs(e * k_wave * phi_amplitude / m_e))
345 v_trap = omega_bounce / k_wave
346
347 print(f"\nTrapping Analysis:")
348 print(f" Wave phase velocity: {v_phase/v_th:.2f} v_th")
349 print(f" Potential amplitude: {phi_amplitude:.2e} V")
350 print(f" Bounce frequency: {omega_bounce:.2e} rad/s")
351 print(f" Trapping width: {v_trap/v_th:.2f} v_th")
352
353 # Create figure
354 fig, axes = plt.subplots(2, 2, figsize=(14, 10))
355
356 # Plot 1: Phase space with separatrix
357 ax = axes[0, 0]
358 cs = ax.contourf(x / lambda_D, v / v_th, f_trap.T, levels=30, cmap='hot')
359 ax.set_xlabel('x [λ_D]', fontsize=12)
360 ax.set_ylabel('v [v_th]', fontsize=12)
361 ax.set_title('Phase Space (Nonlinear Phase)', fontsize=13, fontweight='bold')
362 ax.set_ylim([1, 5])
363 plt.colorbar(cs, ax=ax)
364
365 # Plot 2: Wave potential
366 ax = axes[0, 1]
367 ax.plot(x / lambda_D, phi_trap, 'b-', linewidth=2)
368 ax.set_xlabel('x [λ_D]', fontsize=12)
369 ax.set_ylabel('φ [V]', fontsize=12)
370 ax.set_title('Wave Potential', fontsize=13, fontweight='bold')
371 ax.grid(True, alpha=0.3)
372
373 # Plot 3: Trapped particle orbits (schematic in v-x)
374 ax = axes[1, 0]
375 # Draw approximate separatrix
376 x_sep = np.linspace(0, x[-1], 100)
377 v_sep_upper = v_phase + v_trap * np.cos(k_wave * x_sep)
378 v_sep_lower = v_phase - v_trap * np.cos(k_wave * x_sep)
379
380 ax.plot(x_sep / lambda_D, v_sep_upper / v_th, 'g--', linewidth=2,
381 label='Separatrix (approx)')
382 ax.plot(x_sep / lambda_D, v_sep_lower / v_th, 'g--', linewidth=2)
383 ax.axhline(v_phase / v_th, color='red', linestyle=':', linewidth=2,
384 label=f'v_phase = {v_phase/v_th:.1f} v_th')
385
386 # Overlay phase space contours
387 cs = ax.contour(x / lambda_D, v / v_th, f_trap.T, levels=10,
388 colors='white', linewidths=0.5, alpha=0.5)
389
390 ax.set_xlabel('x [λ_D]', fontsize=12)
391 ax.set_ylabel('v [v_th]', fontsize=12)
392 ax.set_title('Trapping Region (Separatrix)', fontsize=13, fontweight='bold')
393 ax.legend(loc='best', fontsize=10)
394 ax.set_ylim([2, 4])
395
396 # Plot 4: Cut through phase space at x=0
397 ax = axes[1, 1]
398 i_x0 = 0
399 ax.plot(v / v_th, f_trap[i_x0, :] * v_th, 'b-', linewidth=2)
400 ax.axvline(v_phase / v_th, color='red', linestyle='--', linewidth=2,
401 label='Phase velocity')
402 ax.axvspan((v_phase - v_trap) / v_th, (v_phase + v_trap) / v_th,
403 alpha=0.2, color='green', label='Trapping region')
404
405 ax.set_xlabel('v [v_th]', fontsize=12)
406 ax.set_ylabel('f(v) × v_th', fontsize=12)
407 ax.set_title('Velocity Distribution at x=0', fontsize=13, fontweight='bold')
408 ax.legend(loc='best', fontsize=10)
409 ax.grid(True, alpha=0.3)
410 ax.set_xlim([1, 5])
411
412 plt.tight_layout()
413 plt.savefig('bump_on_tail_trapping.png', dpi=300, bbox_inches='tight')
414 plt.show()
415
416
417if __name__ == '__main__':
418 print("\n" + "="*80)
419 print("BUMP-ON-TAIL INSTABILITY: DETAILED ANALYSIS")
420 print("="*80)
421
422 # Run simulation
423 history, solver, omega_pe, lambda_D, v_th = simulate_full_evolution()
424
425 print("\nGenerating analysis plots...")
426 print(" 1. Full evolution (linear → nonlinear → saturation)...")
427 plot_full_evolution(history, solver, omega_pe, lambda_D, v_th)
428
429 print(" 2. Particle trapping analysis...")
430 plot_trapping_analysis(history, omega_pe, lambda_D, v_th)
431
432 print("\nKey Physics:")
433 print(" - Linear phase: Exponential growth due to positive df/dv")
434 print(" - Nonlinear phase: Wave-particle interaction, trapping begins")
435 print(" - Saturation: Plateau formation (df/dv → 0), trapping vortices")
436
437 print("\nDone! Generated 2 figures:")
438 print(" - bump_on_tail_full_evolution.png")
439 print(" - bump_on_tail_trapping.png")