1#!/usr/bin/env python3
2"""
3Electromagnetic Wave Propagation in Unmagnetized Plasma
4
5This script demonstrates EM wave propagation in a plasma, including:
6- Dispersion relation: ω² = ωpe² + k²c²
7- Cutoff at ω = ωpe (evanescent below)
8- Refractive index: n = sqrt(1 - ωpe²/ω²)
9- Wave field visualization entering a plasma density gradient
10
11Key Physics:
12- EM waves cannot propagate below the plasma frequency (cutoff)
13- Evanescent decay with skin depth δ ~ c/ωpe
14- Group velocity vg = c²k/ω < c
15
16Author: Plasma Physics Examples
17License: MIT
18"""
19
20import numpy as np
21import matplotlib.pyplot as plt
22from matplotlib.gridspec import GridSpec
23from matplotlib.animation import FuncAnimation
24
25# Physical constants
26EPS0 = 8.854187817e-12 # F/m
27QE = 1.602176634e-19 # C
28ME = 9.10938356e-31 # kg
29C = 2.99792458e8 # m/s
30
31def compute_plasma_frequency(ne):
32 """
33 Compute electron plasma frequency.
34
35 Parameters:
36 -----------
37 ne : float or array
38 Electron density [m^-3]
39
40 Returns:
41 --------
42 omega_pe : float or array
43 Plasma frequency [rad/s]
44 """
45 omega_pe = np.sqrt(ne * QE**2 / (ME * EPS0))
46 return omega_pe
47
48def em_wave_dispersion(k, omega_pe):
49 """
50 Compute EM wave dispersion in plasma: ω² = ωpe² + k²c².
51
52 Parameters:
53 -----------
54 k : array
55 Wavenumber [rad/m]
56 omega_pe : float
57 Plasma frequency [rad/s]
58
59 Returns:
60 --------
61 omega : array
62 Angular frequency [rad/s]
63 """
64 omega_sq = omega_pe**2 + (k * C)**2
65 return np.sqrt(omega_sq)
66
67def refractive_index(omega, omega_pe):
68 """
69 Compute refractive index: n = sqrt(1 - ωpe²/ω²).
70
71 Parameters:
72 -----------
73 omega : array
74 Angular frequency [rad/s]
75 omega_pe : float
76 Plasma frequency [rad/s]
77
78 Returns:
79 --------
80 n : array
81 Refractive index (complex if ω < ωpe)
82 """
83 ratio = omega_pe**2 / omega**2
84
85 # Real refractive index when ω > ωpe
86 # Imaginary when ω < ωpe (evanescent)
87 n = np.sqrt(1 - ratio + 0j) # Allow complex values
88
89 return n
90
91def skin_depth(omega, omega_pe):
92 """
93 Compute skin depth for evanescent waves (ω < ωpe).
94
95 δ = c / sqrt(ωpe² - ω²)
96
97 Parameters:
98 -----------
99 omega : float
100 Angular frequency [rad/s]
101 omega_pe : float
102 Plasma frequency [rad/s]
103
104 Returns:
105 --------
106 delta : float
107 Skin depth [m]
108 """
109 if omega >= omega_pe:
110 return np.inf
111 else:
112 delta = C / np.sqrt(omega_pe**2 - omega**2)
113 return delta
114
115def simulate_wave_entering_plasma(ne_max, freq, num_cycles=3):
116 """
117 Simulate EM wave entering a plasma density gradient.
118
119 Parameters:
120 -----------
121 ne_max : float
122 Maximum plasma density [m^-3]
123 freq : float
124 Wave frequency [Hz]
125 num_cycles : int
126 Number of wave cycles to simulate
127
128 Returns:
129 --------
130 x, Ex_total : arrays for visualization
131 """
132 omega = 2 * np.pi * freq
133 omega_pe_max = compute_plasma_frequency(ne_max)
134
135 # Spatial grid
136 x = np.linspace(-0.5, 2.0, 2000) # meters
137
138 # Density profile: linear ramp from 0 to ne_max
139 ne = np.zeros_like(x)
140 ramp_start = 0.0
141 ramp_end = 1.0
142 mask = (x >= ramp_start) & (x <= ramp_end)
143 ne[mask] = ne_max * (x[mask] - ramp_start) / (ramp_end - ramp_start)
144 ne[x > ramp_end] = ne_max
145
146 # Plasma frequency profile
147 omega_pe = compute_plasma_frequency(ne)
148
149 # Refractive index profile
150 n = refractive_index(omega, omega_pe)
151
152 # Find cutoff position (where ω = ωpe)
153 cutoff_idx = np.argmin(np.abs(omega - omega_pe[x >= 0]))
154 x_cutoff = x[x >= 0][cutoff_idx] if omega < omega_pe_max else np.inf
155
156 return x, ne, omega_pe, n, x_cutoff
157
158def plot_em_wave_propagation():
159 """
160 Create comprehensive visualization of EM wave propagation in plasma.
161 """
162 # Plasma parameters
163 ne = 1e17 # m^-3 (lower density for radio waves)
164
165 omega_pe = compute_plasma_frequency(ne)
166 f_pe = omega_pe / (2 * np.pi)
167
168 print("=" * 70)
169 print("EM Wave Propagation in Plasma")
170 print("=" * 70)
171 print(f"Electron density: {ne:.2e} m^-3")
172 print(f"Plasma frequency: {f_pe/1e9:.3f} GHz")
173 print(f"Cutoff wavelength: {C/f_pe:.3f} m")
174 print("=" * 70)
175
176 # Create figure
177 fig = plt.figure(figsize=(14, 12))
178 gs = GridSpec(4, 2, figure=fig, hspace=0.35, wspace=0.3)
179
180 # Plot 1: Dispersion relation
181 ax1 = fig.add_subplot(gs[0, :])
182
183 k = np.linspace(0, 3 * omega_pe / C, 1000)
184 omega = em_wave_dispersion(k, omega_pe)
185
186 # Light line in vacuum
187 omega_vacuum = k * C
188
189 ax1.plot(k * C / omega_pe, omega / omega_pe,
190 'b-', linewidth=2, label='Plasma')
191 ax1.plot(k * C / omega_pe, omega_vacuum / omega_pe,
192 'r--', linewidth=2, label='Vacuum (light line)')
193 ax1.axhline(y=1.0, color='green', linestyle=':', linewidth=2,
194 label=r'Cutoff ($\omega = \omega_{pe}$)')
195
196 # Shade forbidden region
197 ax1.fill_between([0, 3], [0, 0], [1, 1], alpha=0.2, color='red',
198 label='Forbidden (evanescent)')
199
200 ax1.set_xlabel(r'$k c / \omega_{pe}$', fontsize=12)
201 ax1.set_ylabel(r'$\omega / \omega_{pe}$', fontsize=12)
202 ax1.set_title('EM Wave Dispersion Relation', fontsize=14, fontweight='bold')
203 ax1.grid(True, alpha=0.3)
204 ax1.legend(fontsize=10, loc='lower right')
205 ax1.set_xlim([0, 3])
206 ax1.set_ylim([0, 3])
207
208 # Plot 2: Refractive index vs frequency
209 ax2 = fig.add_subplot(gs[1, 0])
210
211 omega_array = np.linspace(0.5 * omega_pe, 3 * omega_pe, 1000)
212 n = refractive_index(omega_array, omega_pe)
213
214 ax2.plot(omega_array / omega_pe, np.real(n),
215 'b-', linewidth=2, label='Real part')
216 ax2.plot(omega_array / omega_pe, np.abs(np.imag(n)),
217 'r-', linewidth=2, label='Imaginary part')
218 ax2.axvline(x=1.0, color='green', linestyle=':', linewidth=2,
219 label='Cutoff')
220
221 ax2.set_xlabel(r'$\omega / \omega_{pe}$', fontsize=12)
222 ax2.set_ylabel('Refractive Index n', fontsize=12)
223 ax2.set_title('Refractive Index', fontsize=12, fontweight='bold')
224 ax2.grid(True, alpha=0.3)
225 ax2.legend(fontsize=10)
226 ax2.set_xlim([0.5, 3])
227
228 # Plot 3: Group velocity
229 ax3 = fig.add_subplot(gs[1, 1])
230
231 # For propagating waves (ω > ωpe)
232 omega_prop = omega_array[omega_array > omega_pe]
233 vg = C**2 * np.sqrt(omega_prop**2 - omega_pe**2) / omega_prop
234
235 ax3.plot(omega_prop / omega_pe, vg / C,
236 'b-', linewidth=2)
237 ax3.axvline(x=1.0, color='green', linestyle=':', linewidth=2,
238 label='Cutoff')
239 ax3.axhline(y=1.0, color='r', linestyle='--', linewidth=1,
240 label='Vacuum speed')
241
242 ax3.set_xlabel(r'$\omega / \omega_{pe}$', fontsize=12)
243 ax3.set_ylabel(r'$v_g / c$', fontsize=12)
244 ax3.set_title('Group Velocity', fontsize=12, fontweight='bold')
245 ax3.grid(True, alpha=0.3)
246 ax3.legend(fontsize=10)
247 ax3.set_xlim([1, 3])
248 ax3.set_ylim([0, 1.1])
249
250 # Plot 4: Skin depth vs frequency
251 ax4 = fig.add_subplot(gs[2, 0])
252
253 omega_evan = np.linspace(0.1 * omega_pe, 0.99 * omega_pe, 1000)
254 delta = np.array([skin_depth(w, omega_pe) for w in omega_evan])
255
256 ax4.plot(omega_evan / omega_pe, delta * omega_pe / C,
257 'r-', linewidth=2)
258 ax4.axvline(x=1.0, color='green', linestyle=':', linewidth=2,
259 label='Cutoff')
260
261 ax4.set_xlabel(r'$\omega / \omega_{pe}$', fontsize=12)
262 ax4.set_ylabel(r'Skin Depth ($\omega_{pe} / c$)', fontsize=12)
263 ax4.set_title('Evanescent Wave Skin Depth', fontsize=12, fontweight='bold')
264 ax4.grid(True, alpha=0.3)
265 ax4.legend(fontsize=10)
266 ax4.set_xlim([0.1, 1])
267 ax4.set_yscale('log')
268
269 # Plot 5 & 6: Wave entering plasma gradient
270 # Case 1: Propagating (ω > ωpe)
271 ax5 = fig.add_subplot(gs[2, 1])
272
273 ne_max = 0.5e17 # Lower than ne, so wave can propagate
274 freq_high = 1.5 * f_pe # Above cutoff
275
276 x, ne_profile, omega_pe_profile, n_profile, x_cutoff = \
277 simulate_wave_entering_plasma(ne_max, freq_high)
278
279 ax5_twin = ax5.twinx()
280 ax5.plot(x, ne_profile / 1e17, 'g-', linewidth=2, label='Density')
281 ax5.set_xlabel('Position (m)', fontsize=10)
282 ax5.set_ylabel(r'$n_e$ ($10^{17}$ m$^{-3}$)', fontsize=10, color='g')
283 ax5.tick_params(axis='y', labelcolor='g')
284
285 ax5_twin.plot(x, np.real(n_profile), 'b-', linewidth=2, label='Re(n)')
286 ax5_twin.set_ylabel('Refractive Index', fontsize=10, color='b')
287 ax5_twin.tick_params(axis='y', labelcolor='b')
288 ax5_twin.axhline(y=0, color='r', linestyle='--', linewidth=1)
289
290 ax5.set_title(f'Propagating Wave (f = {freq_high/1e9:.2f} GHz > fpe)',
291 fontsize=11, fontweight='bold')
292 ax5.grid(True, alpha=0.3)
293 ax5.set_xlim([-0.5, 2])
294
295 # Case 2: Evanescent (ω < ωpe)
296 ax6 = fig.add_subplot(gs[3, :])
297
298 ne_max = 2e17 # Higher than ne, wave will be cutoff
299 freq_low = 0.7 * f_pe # Below cutoff
300
301 x, ne_profile, omega_pe_profile, n_profile, x_cutoff = \
302 simulate_wave_entering_plasma(ne_max, freq_low)
303
304 ax6_twin = ax6.twinx()
305 ax6.plot(x, ne_profile / 1e17, 'g-', linewidth=2, label='Density')
306 ax6.axvline(x=x_cutoff, color='purple', linestyle=':', linewidth=2,
307 label=f'Cutoff at x = {x_cutoff:.3f} m')
308 ax6.set_xlabel('Position (m)', fontsize=10)
309 ax6.set_ylabel(r'$n_e$ ($10^{17}$ m$^{-3}$)', fontsize=10, color='g')
310 ax6.tick_params(axis='y', labelcolor='g')
311
312 ax6_twin.plot(x, np.real(n_profile), 'b-', linewidth=2, label='Re(n)')
313 ax6_twin.plot(x, np.abs(np.imag(n_profile)), 'r--', linewidth=2, label='|Im(n)|')
314 ax6_twin.set_ylabel('Refractive Index', fontsize=10, color='b')
315 ax6_twin.tick_params(axis='y', labelcolor='b')
316 ax6_twin.axhline(y=0, color='gray', linestyle='-', linewidth=0.5)
317
318 ax6.set_title(f'Evanescent Wave (f = {freq_low/1e9:.2f} GHz < fpe)',
319 fontsize=11, fontweight='bold')
320 ax6.grid(True, alpha=0.3)
321 ax6.legend(loc='upper left', fontsize=9)
322 ax6_twin.legend(loc='upper right', fontsize=9)
323 ax6.set_xlim([-0.5, 2])
324
325 plt.suptitle('Electromagnetic Wave Propagation in Unmagnetized Plasma',
326 fontsize=16, fontweight='bold', y=0.997)
327
328 plt.savefig('em_wave_plasma.png', dpi=150, bbox_inches='tight')
329 print("\nPlot saved as 'em_wave_plasma.png'")
330
331 plt.show()
332
333if __name__ == "__main__":
334 plot_em_wave_propagation()