1#!/usr/bin/env python3
2"""
3Langmuir Wave Dispersion Relation
4
5This script computes and visualizes the electrostatic Langmuir wave dispersion
6relation in both cold and warm (kinetic) plasmas. It demonstrates the Bohm-Gross
7dispersion relation and shows where Landau damping becomes important.
8
9Key Physics:
10- Cold plasma: ω = ωpe (no dependence on k)
11- Warm plasma (Bohm-Gross): ω² = ωpe² + 3k²vth²
12- Landau damping becomes significant when k*λD ~ 1
13
14Author: Plasma Physics Examples
15License: MIT
16"""
17
18import numpy as np
19import matplotlib.pyplot as plt
20from matplotlib.gridspec import GridSpec
21
22# Physical constants
23EPS0 = 8.854187817e-12 # F/m
24QE = 1.602176634e-19 # C
25ME = 9.10938356e-31 # kg
26KB = 1.380649e-23 # J/K
27
28def compute_plasma_parameters(ne, Te):
29 """
30 Compute fundamental plasma parameters.
31
32 Parameters:
33 -----------
34 ne : float
35 Electron density [m^-3]
36 Te : float
37 Electron temperature [eV]
38
39 Returns:
40 --------
41 dict : Dictionary containing plasma parameters
42 """
43 # Convert temperature to SI units
44 Te_joule = Te * QE
45
46 # Electron plasma frequency
47 omega_pe = np.sqrt(ne * QE**2 / (ME * EPS0))
48 f_pe = omega_pe / (2 * np.pi)
49
50 # Electron thermal velocity
51 vth = np.sqrt(2 * Te_joule / ME)
52
53 # Debye length
54 lambda_D = np.sqrt(EPS0 * Te_joule / (ne * QE**2))
55
56 return {
57 'omega_pe': omega_pe,
58 'f_pe': f_pe,
59 'vth': vth,
60 'lambda_D': lambda_D,
61 'ne': ne,
62 'Te': Te
63 }
64
65def bohm_gross_dispersion(k, params):
66 """
67 Compute Bohm-Gross dispersion relation for Langmuir waves.
68
69 ω² = ωpe² + 3k²vth²
70
71 Parameters:
72 -----------
73 k : array
74 Wavenumber [rad/m]
75 params : dict
76 Plasma parameters
77
78 Returns:
79 --------
80 omega : array
81 Angular frequency [rad/s]
82 """
83 omega_pe = params['omega_pe']
84 vth = params['vth']
85
86 omega_sq = omega_pe**2 + 3 * k**2 * vth**2
87 return np.sqrt(omega_sq)
88
89def cold_plasma_dispersion(k, params):
90 """
91 Cold plasma dispersion: ω = ωpe (constant).
92
93 Parameters:
94 -----------
95 k : array
96 Wavenumber [rad/m]
97 params : dict
98 Plasma parameters
99
100 Returns:
101 --------
102 omega : array
103 Angular frequency [rad/s]
104 """
105 omega_pe = params['omega_pe']
106 return omega_pe * np.ones_like(k)
107
108def phase_velocity(omega, k):
109 """
110 Compute phase velocity vph = ω/k.
111
112 Parameters:
113 -----------
114 omega : array
115 Angular frequency [rad/s]
116 k : array
117 Wavenumber [rad/m]
118
119 Returns:
120 --------
121 vph : array
122 Phase velocity [m/s]
123 """
124 return omega / k
125
126def group_velocity(omega, k):
127 """
128 Compute group velocity vg = dω/dk numerically.
129
130 Parameters:
131 -----------
132 omega : array
133 Angular frequency [rad/s]
134 k : array
135 Wavenumber [rad/m]
136
137 Returns:
138 --------
139 vg : array
140 Group velocity [m/s]
141 """
142 # Use central differences for interior points
143 vg = np.gradient(omega, k)
144 return vg
145
146def plot_langmuir_dispersion():
147 """
148 Create comprehensive visualization of Langmuir wave dispersion.
149 """
150 # Define plasma parameters for a typical laboratory plasma
151 ne = 1e18 # m^-3 (10^12 cm^-3)
152 Te = 10.0 # eV
153
154 params = compute_plasma_parameters(ne, Te)
155
156 # Print plasma parameters
157 print("=" * 60)
158 print("Plasma Parameters")
159 print("=" * 60)
160 print(f"Electron density: {params['ne']:.2e} m^-3")
161 print(f"Electron temperature: {params['Te']:.1f} eV")
162 print(f"Electron plasma frequency: {params['f_pe']/1e9:.3f} GHz")
163 print(f"Thermal velocity: {params['vth']/1e6:.3f} × 10^6 m/s")
164 print(f"Debye length: {params['lambda_D']*1e6:.3f} μm")
165 print("=" * 60)
166
167 # Create wavenumber array
168 # Range from small k to k ~ 10/λD
169 k_min = 1e3 # rad/m
170 k_max = 10 / params['lambda_D']
171 k = np.linspace(k_min, k_max, 1000)
172
173 # Compute dispersion relations
174 omega_warm = bohm_gross_dispersion(k, params)
175 omega_cold = cold_plasma_dispersion(k, params)
176
177 # Compute velocities
178 vph_warm = phase_velocity(omega_warm, k)
179 vg_warm = group_velocity(omega_warm, k)
180
181 # Identify Landau damping region (k*λD ~ 0.5 to 1.5)
182 k_lambda_D = k * params['lambda_D']
183 landau_mask = (k_lambda_D >= 0.5) & (k_lambda_D <= 1.5)
184
185 # Create figure with subplots
186 fig = plt.figure(figsize=(14, 10))
187 gs = GridSpec(3, 2, figure=fig, hspace=0.3, wspace=0.3)
188
189 # Plot 1: Dispersion relation ω(k)
190 ax1 = fig.add_subplot(gs[0, :])
191 ax1.plot(k * params['lambda_D'], omega_warm / params['omega_pe'],
192 'b-', linewidth=2, label='Warm plasma (Bohm-Gross)')
193 ax1.plot(k * params['lambda_D'], omega_cold / params['omega_pe'],
194 'r--', linewidth=2, label='Cold plasma')
195
196 # Highlight Landau damping region
197 ax1.axvspan(0.5, 1.5, alpha=0.2, color='yellow',
198 label='Strong Landau damping (k·λD ~ 1)')
199
200 ax1.set_xlabel(r'$k \lambda_D$', fontsize=12)
201 ax1.set_ylabel(r'$\omega / \omega_{pe}$', fontsize=12)
202 ax1.set_title('Langmuir Wave Dispersion Relation', fontsize=14, fontweight='bold')
203 ax1.grid(True, alpha=0.3)
204 ax1.legend(fontsize=10)
205 ax1.set_xlim([0, 10])
206
207 # Plot 2: Phase velocity
208 ax2 = fig.add_subplot(gs[1, 0])
209 ax2.plot(k * params['lambda_D'], vph_warm / params['vth'],
210 'b-', linewidth=2, label='Warm plasma')
211 ax2.axhline(y=0, color='r', linestyle='--', linewidth=2, label='Cold plasma (∞)')
212
213 # Highlight Landau damping region
214 ax2.axvspan(0.5, 1.5, alpha=0.2, color='yellow')
215
216 ax2.set_xlabel(r'$k \lambda_D$', fontsize=12)
217 ax2.set_ylabel(r'$v_{ph} / v_{th}$', fontsize=12)
218 ax2.set_title('Phase Velocity', fontsize=12, fontweight='bold')
219 ax2.grid(True, alpha=0.3)
220 ax2.legend(fontsize=10)
221 ax2.set_xlim([0, 10])
222 ax2.set_ylim([0, 20])
223
224 # Plot 3: Group velocity
225 ax3 = fig.add_subplot(gs[1, 1])
226 ax3.plot(k * params['lambda_D'], vg_warm / params['vth'],
227 'b-', linewidth=2, label='Warm plasma')
228 ax3.axhline(y=0, color='r', linestyle='--', linewidth=2, label='Cold plasma')
229
230 # Highlight Landau damping region
231 ax3.axvspan(0.5, 1.5, alpha=0.2, color='yellow')
232
233 ax3.set_xlabel(r'$k \lambda_D$', fontsize=12)
234 ax3.set_ylabel(r'$v_g / v_{th}$', fontsize=12)
235 ax3.set_title('Group Velocity', fontsize=12, fontweight='bold')
236 ax3.grid(True, alpha=0.3)
237 ax3.legend(fontsize=10)
238 ax3.set_xlim([0, 10])
239
240 # Plot 4: Comparison in frequency space
241 ax4 = fig.add_subplot(gs[2, 0])
242 f_warm = omega_warm / (2 * np.pi * params['f_pe'])
243 f_cold = omega_cold / (2 * np.pi * params['f_pe'])
244
245 ax4.plot(k * params['lambda_D'], f_warm,
246 'b-', linewidth=2, label='Warm plasma')
247 ax4.plot(k * params['lambda_D'], f_cold,
248 'r--', linewidth=2, label='Cold plasma')
249
250 # Highlight Landau damping region
251 ax4.axvspan(0.5, 1.5, alpha=0.2, color='yellow')
252
253 ax4.set_xlabel(r'$k \lambda_D$', fontsize=12)
254 ax4.set_ylabel(r'$f / f_{pe}$', fontsize=12)
255 ax4.set_title('Frequency Normalization', fontsize=12, fontweight='bold')
256 ax4.grid(True, alpha=0.3)
257 ax4.legend(fontsize=10)
258 ax4.set_xlim([0, 10])
259
260 # Plot 5: Relative correction
261 ax5 = fig.add_subplot(gs[2, 1])
262 correction = (omega_warm - omega_cold) / omega_cold * 100
263
264 ax5.plot(k * params['lambda_D'], correction,
265 'g-', linewidth=2)
266
267 # Highlight Landau damping region
268 ax5.axvspan(0.5, 1.5, alpha=0.2, color='yellow')
269
270 ax5.set_xlabel(r'$k \lambda_D$', fontsize=12)
271 ax5.set_ylabel('Correction (%)', fontsize=12)
272 ax5.set_title('Kinetic Correction to Cold Plasma Frequency',
273 fontsize=12, fontweight='bold')
274 ax5.grid(True, alpha=0.3)
275 ax5.set_xlim([0, 10])
276
277 plt.suptitle('Langmuir Wave Characteristics in Warm vs Cold Plasma',
278 fontsize=16, fontweight='bold', y=0.995)
279
280 plt.savefig('langmuir_wave_dispersion.png', dpi=150, bbox_inches='tight')
281 print("\nPlot saved as 'langmuir_wave_dispersion.png'")
282
283 plt.show()
284
285if __name__ == "__main__":
286 plot_langmuir_dispersion()