1#!/usr/bin/env python3
2"""
3Two-Fluid vs MHD Wave Dispersion
4
5This script compares single-fluid MHD wave dispersion with two-fluid
6corrections for parallel propagation in a magnetized plasma.
7
8Key Physics:
9- MHD AlfvĆ©n wave: Ļ = kĀ·vA (linear dispersion)
10- Two-fluid corrections at kĀ·di ~ 1 (ion skin depth)
11- Ion cyclotron wave and whistler branches
12- Kinetic AlfvƩn wave at finite k_perp
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
26MP = 1.672621898e-27 # kg
27C = 2.99792458e8 # m/s
28MU0 = 4 * np.pi * 1e-7 # H/m
29
30def compute_plasma_parameters(ne, B0, mi=MP):
31 """
32 Compute characteristic plasma parameters.
33
34 Parameters:
35 -----------
36 ne : float
37 Electron density [m^-3]
38 B0 : float
39 Magnetic field [T]
40 mi : float
41 Ion mass [kg]
42
43 Returns:
44 --------
45 dict : Plasma parameters
46 """
47 # AlfvƩn speed
48 vA = B0 / np.sqrt(MU0 * ne * mi)
49
50 # Ion skin depth
51 di = C / np.sqrt(ne * QE**2 / (mi * EPS0))
52
53 # Cyclotron frequencies
54 omega_ci = QE * B0 / mi
55 omega_ce = QE * B0 / ME
56
57 # Plasma frequencies
58 omega_pi = np.sqrt(ne * QE**2 / (mi * EPS0))
59 omega_pe = np.sqrt(ne * QE**2 / (ME * EPS0))
60
61 # Ion Larmor radius (thermal)
62 # Assume Ti = 1 keV for estimate
63 Ti = 1000 * QE # J
64 vti = np.sqrt(2 * Ti / mi)
65 rho_i = vti / omega_ci
66
67 return {
68 'vA': vA,
69 'di': di,
70 'rho_i': rho_i,
71 'omega_ci': omega_ci,
72 'omega_ce': omega_ce,
73 'omega_pi': omega_pi,
74 'omega_pe': omega_pe,
75 'ne': ne,
76 'B0': B0
77 }
78
79def mhd_alfven_dispersion(k, params):
80 """
81 MHD AlfvĆ©n wave: Ļ = kĀ·vA.
82
83 Parameters:
84 -----------
85 k : array
86 Wavenumber [rad/m]
87 params : dict
88 Plasma parameters
89
90 Returns:
91 --------
92 omega : array
93 Angular frequency [rad/s]
94 """
95 return k * params['vA']
96
97def two_fluid_parallel_dispersion(k, params):
98 """
99 Two-fluid dispersion for parallel propagation.
100
101 Includes ion cyclotron and whistler branches.
102
103 Approximate dispersion:
104 ĻĀ² = kĀ²vAĀ² * (1 + kĀ²diĀ²) / (1 + kĀ²diĀ² + vAĀ²/cĀ²)
105
106 Parameters:
107 -----------
108 k : array
109 Wavenumber [rad/m]
110 params : dict
111 Plasma parameters
112
113 Returns:
114 --------
115 omega : array
116 Angular frequency [rad/s]
117 """
118 vA = params['vA']
119 di = params['di']
120 omega_ci = params['omega_ci']
121
122 # Two-fluid correction
123 k_di = k * di
124
125 # Low frequency branch (ion cyclotron)
126 # Ļ ā kĀ·vA / sqrt(1 + kĀ²diĀ²)
127 omega_ic = k * vA / np.sqrt(1 + k_di**2)
128
129 return omega_ic
130
131def whistler_branch(k, params):
132 """
133 High-frequency whistler branch in two-fluid theory.
134
135 Ļ ā (kĀ²vAĀ²)/(ĻciĀ·diĀ·k) for kĀ·di >> 1
136
137 Parameters:
138 -----------
139 k : array
140 Wavenumber [rad/m]
141 params : dict
142 Plasma parameters
143
144 Returns:
145 --------
146 omega : array
147 Angular frequency [rad/s]
148 """
149 vA = params['vA']
150 di = params['di']
151 omega_ci = params['omega_ci']
152 omega_ce = params['omega_ce']
153
154 # Whistler approximation: Ļ ā kĀ²vAĀ²/(Ļci) for kĀ·di >> 1
155 # More accurate: Ļ ā kĀ²cĀ²/(Ļpe) * (Ļce/Ļ)
156 # Simple form:
157 omega_w = (k * vA)**2 / (omega_ci + 1e-10) # Avoid division by zero
158
159 # Cap at electron cyclotron frequency
160 omega_w = np.minimum(omega_w, 0.5 * omega_ce)
161
162 return omega_w
163
164def kinetic_alfven_dispersion(k_perp, k_parallel, params):
165 """
166 Kinetic AlfvƩn wave dispersion with finite k_perp.
167
168 Ļ = k_ā„Ā·vA * sqrt(1 + k_perpĀ²Ā·ĻiĀ²)
169
170 Parameters:
171 -----------
172 k_perp : array
173 Perpendicular wavenumber [rad/m]
174 k_parallel : float
175 Parallel wavenumber [rad/m]
176 params : dict
177 Plasma parameters
178
179 Returns:
180 --------
181 omega : array
182 Angular frequency [rad/s]
183 """
184 vA = params['vA']
185 rho_i = params['rho_i']
186
187 # KAW dispersion
188 omega = k_parallel * vA * np.sqrt(1 + (k_perp * rho_i)**2)
189
190 return omega
191
192def plot_two_fluid_waves():
193 """
194 Create comprehensive comparison of MHD vs two-fluid waves.
195 """
196 # Typical tokamak parameters
197 ne = 1e19 # m^-3
198 B0 = 2.0 # T
199 mi = 2 * MP # Deuterium
200
201 params = compute_plasma_parameters(ne, B0, mi)
202
203 print("=" * 70)
204 print("Two-Fluid vs MHD Wave Comparison")
205 print("=" * 70)
206 print(f"Electron density: {params['ne']:.2e} m^-3")
207 print(f"Magnetic field: {params['B0']:.2f} T")
208 print(f"AlfvĆ©n speed: {params['vA']/1e6:.3f} Ć 10^6 m/s")
209 print(f"Ion skin depth: {params['di']*100:.2f} cm")
210 print(f"Ion Larmor radius (1 keV): {params['rho_i']*1000:.2f} mm")
211 print(f"Ion cyclotron frequency: {params['omega_ci']/(2*np.pi)/1e6:.2f} MHz")
212 print("=" * 70)
213
214 # Create figure
215 fig = plt.figure(figsize=(14, 12))
216 gs = GridSpec(3, 2, figure=fig, hspace=0.35, wspace=0.3)
217
218 # Wavenumber array
219 k_min = 1e0
220 k_max = 100 / params['di']
221 k = np.logspace(np.log10(k_min), np.log10(k_max), 1000)
222
223 # Compute dispersions
224 omega_mhd = mhd_alfven_dispersion(k, params)
225 omega_2f = two_fluid_parallel_dispersion(k, params)
226 omega_whistler = whistler_branch(k, params)
227
228 # Normalize
229 k_di = k * params['di']
230 omega_norm_mhd = omega_mhd / params['omega_ci']
231 omega_norm_2f = omega_2f / params['omega_ci']
232 omega_norm_w = omega_whistler / params['omega_ci']
233
234 # Plot 1: Dispersion relation (full range)
235 ax1 = fig.add_subplot(gs[0, :])
236
237 ax1.loglog(k_di, omega_norm_mhd, 'r--', linewidth=2.5, label='MHD AlfvƩn')
238 ax1.loglog(k_di, omega_norm_2f, 'b-', linewidth=2.5, label='Two-fluid (IC branch)')
239 ax1.loglog(k_di, omega_norm_w, 'g-', linewidth=2.5, label='Whistler branch')
240
241 # Mark transition region
242 ax1.axvline(x=1, color='purple', linestyle=':', linewidth=2,
243 label=r'$k \cdot d_i = 1$')
244 ax1.axhline(y=1, color='orange', linestyle=':', linewidth=2,
245 label=r'$\omega = \Omega_{ci}$')
246
247 ax1.set_xlabel(r'$k \cdot d_i$', fontsize=12)
248 ax1.set_ylabel(r'$\omega / \Omega_{ci}$', fontsize=12)
249 ax1.set_title('Parallel Propagation: MHD vs Two-Fluid',
250 fontsize=14, fontweight='bold')
251 ax1.grid(True, alpha=0.3, which='both')
252 ax1.legend(fontsize=10, loc='lower right')
253 ax1.set_xlim([1e-2, 1e2])
254 ax1.set_ylim([1e-2, 1e3])
255
256 # Plot 2: Phase velocity
257 ax2 = fig.add_subplot(gs[1, 0])
258
259 vph_mhd = omega_mhd / k
260 vph_2f = omega_2f / k
261
262 ax2.semilogx(k_di, vph_mhd / params['vA'], 'r--', linewidth=2, label='MHD')
263 ax2.semilogx(k_di, vph_2f / params['vA'], 'b-', linewidth=2, label='Two-fluid')
264
265 ax2.axvline(x=1, color='purple', linestyle=':', linewidth=2)
266 ax2.axhline(y=1, color='gray', linestyle='--', linewidth=1)
267
268 ax2.set_xlabel(r'$k \cdot d_i$', fontsize=12)
269 ax2.set_ylabel(r'$v_{ph} / v_A$', fontsize=12)
270 ax2.set_title('Phase Velocity', fontsize=12, fontweight='bold')
271 ax2.grid(True, alpha=0.3)
272 ax2.legend(fontsize=10)
273 ax2.set_xlim([1e-2, 1e2])
274
275 # Plot 3: Group velocity
276 ax3 = fig.add_subplot(gs[1, 1])
277
278 vg_2f = np.gradient(omega_2f, k)
279
280 ax3.semilogx(k_di, params['vA'] * np.ones_like(k) / params['vA'],
281 'r--', linewidth=2, label='MHD')
282 ax3.semilogx(k_di, vg_2f / params['vA'], 'b-', linewidth=2, label='Two-fluid')
283
284 ax3.axvline(x=1, color='purple', linestyle=':', linewidth=2)
285 ax3.axhline(y=1, color='gray', linestyle='--', linewidth=1)
286
287 ax3.set_xlabel(r'$k \cdot d_i$', fontsize=12)
288 ax3.set_ylabel(r'$v_g / v_A$', fontsize=12)
289 ax3.set_title('Group Velocity', fontsize=12, fontweight='bold')
290 ax3.grid(True, alpha=0.3)
291 ax3.legend(fontsize=10)
292 ax3.set_xlim([1e-2, 1e2])
293
294 # Plot 4: Dispersion relation (linear scale, low k)
295 ax4 = fig.add_subplot(gs[2, 0])
296
297 k_low = np.linspace(0, 5 / params['di'], 500)
298 omega_mhd_low = mhd_alfven_dispersion(k_low, params)
299 omega_2f_low = two_fluid_parallel_dispersion(k_low, params)
300
301 ax4.plot(k_low * params['di'], omega_mhd_low / params['omega_ci'],
302 'r--', linewidth=2, label='MHD')
303 ax4.plot(k_low * params['di'], omega_2f_low / params['omega_ci'],
304 'b-', linewidth=2, label='Two-fluid')
305
306 ax4.axvline(x=1, color='purple', linestyle=':', linewidth=2,
307 label=r'$k \cdot d_i = 1$')
308
309 ax4.set_xlabel(r'$k \cdot d_i$', fontsize=12)
310 ax4.set_ylabel(r'$\omega / \Omega_{ci}$', fontsize=12)
311 ax4.set_title('Low-k Regime (Linear Scale)', fontsize=12, fontweight='bold')
312 ax4.grid(True, alpha=0.3)
313 ax4.legend(fontsize=10)
314
315 # Plot 5: Kinetic AlfvƩn wave (perpendicular k)
316 ax5 = fig.add_subplot(gs[2, 1])
317
318 k_parallel_fixed = 1.0 / params['di'] # Fixed k_ā„
319 k_perp = np.linspace(0, 10 / params['rho_i'], 500)
320
321 omega_kaw = kinetic_alfven_dispersion(k_perp, k_parallel_fixed, params)
322 omega_mhd_perp = k_parallel_fixed * params['vA'] * np.ones_like(k_perp)
323
324 ax5.plot(k_perp * params['rho_i'], omega_kaw / params['omega_ci'],
325 'b-', linewidth=2, label='KAW')
326 ax5.plot(k_perp * params['rho_i'], omega_mhd_perp / params['omega_ci'],
327 'r--', linewidth=2, label='MHD (k_perp = 0)')
328
329 ax5.axvline(x=1, color='green', linestyle=':', linewidth=2,
330 label=r'$k_\perp \cdot \rho_i = 1$')
331
332 ax5.set_xlabel(r'$k_\perp \cdot \rho_i$', fontsize=12)
333 ax5.set_ylabel(r'$\omega / \Omega_{ci}$', fontsize=12)
334 ax5.set_title(f'Kinetic AlfvĆ©n Wave (kā„Ā·di = {k_parallel_fixed * params["di"]:.1f})',
335 fontsize=12, fontweight='bold')
336 ax5.grid(True, alpha=0.3)
337 ax5.legend(fontsize=10)
338
339 plt.suptitle('Two-Fluid Corrections to MHD Wave Dispersion',
340 fontsize=16, fontweight='bold', y=0.997)
341
342 plt.savefig('two_fluid_waves.png', dpi=150, bbox_inches='tight')
343 print("\nPlot saved as 'two_fluid_waves.png'")
344 print("\nKey findings:")
345 print(f" - MHD valid for kĀ·di << 1")
346 print(f" - Two-fluid effects important at kĀ·di ~ 1")
347 print(f" - KAW effects at k_perpĀ·Ļi ~ 1")
348
349 plt.show()
350
351if __name__ == "__main__":
352 plot_two_fluid_waves()