1#!/usr/bin/env python3
2"""
3Collision Frequencies and Transport in Plasmas
4
5This script calculates and visualizes collision frequencies, Coulomb logarithm,
6Spitzer resistivity, and mean free paths for various plasma conditions.
7
8Author: Plasma Physics Examples
9License: MIT
10"""
11
12import numpy as np
13import matplotlib.pyplot as plt
14from scipy.constants import e, m_e, m_p, epsilon_0, k, c
15
16
17def coulomb_logarithm(n_e, T_e, Z=1):
18 """
19 Calculate the Coulomb logarithm.
20
21 Parameters:
22 -----------
23 n_e : float or array
24 Electron density [m^-3]
25 T_e : float or array
26 Electron temperature [eV]
27 Z : int
28 Ion charge number
29
30 Returns:
31 --------
32 float or array : Coulomb logarithm (dimensionless)
33 """
34 # Two regimes based on temperature
35 # High temperature (T_e > 10 eV): classical formula
36 # Low temperature: quantum effects important
37
38 ln_Lambda = np.where(
39 T_e > 10,
40 24 - np.log(np.sqrt(n_e * 1e-6) / T_e), # High T
41 23 - np.log(np.sqrt(n_e * 1e-6) * Z / T_e**1.5) # Low T
42 )
43
44 # Enforce lower bound
45 ln_Lambda = np.maximum(ln_Lambda, 2.0)
46
47 return ln_Lambda
48
49
50def electron_electron_collision_freq(n_e, T_e):
51 """
52 Electron-electron collision frequency.
53
54 Parameters:
55 -----------
56 n_e : float or array
57 Electron density [m^-3]
58 T_e : float or array
59 Electron temperature [eV]
60
61 Returns:
62 --------
63 float or array : ν_ee [Hz]
64 """
65 T_J = T_e * e
66 ln_Lambda = coulomb_logarithm(n_e, T_e)
67
68 # Formula from NRL Plasma Formulary
69 nu_ee = (n_e * e**4 * ln_Lambda) / (
70 12 * np.pi**1.5 * epsilon_0**2 * m_e**0.5 * T_J**1.5
71 )
72
73 return nu_ee
74
75
76def electron_ion_collision_freq(n_e, T_e, Z=1):
77 """
78 Electron-ion collision frequency (momentum transfer).
79
80 Parameters:
81 -----------
82 n_e : float or array
83 Electron density [m^-3]
84 T_e : float or array
85 Electron temperature [eV]
86 Z : int
87 Ion charge number
88
89 Returns:
90 --------
91 float or array : ν_ei [Hz]
92 """
93 T_J = T_e * e
94 ln_Lambda = coulomb_logarithm(n_e, T_e, Z)
95
96 # Spitzer formula
97 nu_ei = (Z * n_e * e**4 * ln_Lambda) / (
98 12 * np.pi**1.5 * epsilon_0**2 * m_e**0.5 * T_J**1.5
99 )
100
101 return nu_ei
102
103
104def ion_ion_collision_freq(n_i, T_i, A=1, Z=1):
105 """
106 Ion-ion collision frequency.
107
108 Parameters:
109 -----------
110 n_i : float or array
111 Ion density [m^-3]
112 T_i : float or array
113 Ion temperature [eV]
114 A : float
115 Ion mass number
116 Z : int
117 Ion charge number
118
119 Returns:
120 --------
121 float or array : ν_ii [Hz]
122 """
123 T_J = T_i * e
124 m_i = A * m_p
125 ln_Lambda = coulomb_logarithm(n_i, T_i, Z)
126
127 nu_ii = (Z**4 * n_i * e**4 * ln_Lambda) / (
128 12 * np.pi**1.5 * epsilon_0**2 * m_i**0.5 * T_J**1.5
129 )
130
131 return nu_ii
132
133
134def spitzer_resistivity(n_e, T_e, Z=1):
135 """
136 Spitzer resistivity (classical).
137
138 Parameters:
139 -----------
140 n_e : float or array
141 Electron density [m^-3]
142 T_e : float or array
143 Electron temperature [eV]
144 Z : int
145 Ion charge number
146
147 Returns:
148 --------
149 float or array : η [Ω·m]
150 """
151 T_J = T_e * e
152 ln_Lambda = coulomb_logarithm(n_e, T_e, Z)
153
154 # Spitzer formula with numerical factor
155 eta = (Z * m_e * e**2 * ln_Lambda) / (
156 12 * np.pi**1.5 * epsilon_0**2 * T_J**1.5
157 )
158
159 return eta
160
161
162def mean_free_path(n_e, T_e):
163 """
164 Electron mean free path.
165
166 Parameters:
167 -----------
168 n_e : float or array
169 Electron density [m^-3]
170 T_e : float or array
171 Electron temperature [eV]
172
173 Returns:
174 --------
175 float or array : λ_mfp [m]
176 """
177 T_J = T_e * e
178 v_th = np.sqrt(2 * T_J / m_e)
179 nu_ei = electron_ion_collision_freq(n_e, T_e)
180
181 return v_th / nu_ei
182
183
184def plot_collision_frequencies():
185 """Plot collision frequencies as function of temperature."""
186
187 n_e = 1e20 # m^-3 (tokamak)
188 T_range = np.logspace(0, 5, 200) # 1 eV to 100 keV
189
190 # Calculate frequencies
191 nu_ee = electron_electron_collision_freq(n_e, T_range)
192 nu_ei = electron_ion_collision_freq(n_e, T_range)
193 nu_ii = ion_ion_collision_freq(n_e, T_range)
194
195 # Plasma frequency for reference
196 omega_pe = np.sqrt(n_e * e**2 / (epsilon_0 * m_e))
197 f_pe = omega_pe / (2 * np.pi)
198
199 # Gyrofrequency (assume B = 5 T)
200 B = 5.0 # Tesla
201 omega_ce = e * B / m_e
202 f_ce = omega_ce / (2 * np.pi)
203
204 # Create figure
205 fig, ax = plt.subplots(figsize=(10, 8))
206
207 ax.loglog(T_range, nu_ee, 'b-', linewidth=2.5, label='ν_ee (e-e collision)')
208 ax.loglog(T_range, nu_ei, 'r-', linewidth=2.5, label='ν_ei (e-i collision)')
209 ax.loglog(T_range, nu_ii, 'g-', linewidth=2.5, label='ν_ii (i-i collision)')
210
211 # Reference frequencies
212 ax.axhline(f_pe, color='purple', linestyle='--', linewidth=2,
213 label=f'f_pe = {f_pe:.2e} Hz')
214 ax.axhline(f_ce, color='orange', linestyle='--', linewidth=2,
215 label=f'f_ce = {f_ce:.2e} Hz (B=5T)')
216
217 ax.set_xlabel('Electron Temperature [eV]', fontsize=13)
218 ax.set_ylabel('Collision Frequency [Hz]', fontsize=13)
219 ax.set_title(f'Collision Frequencies vs Temperature\n(n_e = {n_e:.1e} m⁻³)',
220 fontsize=14, fontweight='bold')
221 ax.legend(loc='best', fontsize=11)
222 ax.grid(True, alpha=0.3, which='both')
223
224 # Add regime annotations
225 ax.text(2, 1e6, 'Collisional\nRegime', fontsize=11, color='red',
226 ha='left', va='center', bbox=dict(boxstyle='round', facecolor='white', alpha=0.7))
227 ax.text(1e4, 1e6, 'Collisionless\nRegime', fontsize=11, color='blue',
228 ha='left', va='center', bbox=dict(boxstyle='round', facecolor='white', alpha=0.7))
229
230 plt.tight_layout()
231 plt.savefig('collision_frequencies.png', dpi=300, bbox_inches='tight')
232 plt.show()
233
234
235def plot_coulomb_logarithm():
236 """Plot Coulomb logarithm for various conditions."""
237
238 # Temperature range
239 T_range = np.logspace(-1, 5, 100) # 0.1 eV to 100 keV
240
241 # Different densities
242 densities = [1e14, 1e16, 1e18, 1e20, 1e22] # m^-3
243
244 fig, ax = plt.subplots(figsize=(10, 7))
245
246 colors = plt.cm.viridis(np.linspace(0.1, 0.9, len(densities)))
247
248 for n_e, color in zip(densities, colors):
249 ln_Lambda = coulomb_logarithm(n_e * np.ones_like(T_range), T_range)
250 ax.semilogx(T_range, ln_Lambda, linewidth=2.5, color=color,
251 label=f'n_e = {n_e:.1e} m⁻³')
252
253 ax.set_xlabel('Temperature [eV]', fontsize=13)
254 ax.set_ylabel('Coulomb Logarithm ln Λ', fontsize=13)
255 ax.set_title('Coulomb Logarithm vs Temperature and Density',
256 fontsize=14, fontweight='bold')
257 ax.legend(loc='best', fontsize=10)
258 ax.grid(True, alpha=0.3)
259 ax.set_ylim([0, 30])
260
261 plt.tight_layout()
262 plt.savefig('coulomb_logarithm.png', dpi=300, bbox_inches='tight')
263 plt.show()
264
265
266def plot_spitzer_resistivity():
267 """Plot Spitzer resistivity vs temperature."""
268
269 # Temperature range
270 T_range = np.logspace(0, 5, 200) # 1 eV to 100 keV
271
272 # Different densities
273 densities = [1e18, 1e19, 1e20, 1e21] # m^-3
274
275 fig, ax = plt.subplots(figsize=(10, 7))
276
277 colors = plt.cm.plasma(np.linspace(0.1, 0.9, len(densities)))
278
279 for n_e, color in zip(densities, colors):
280 eta = spitzer_resistivity(n_e * np.ones_like(T_range), T_range)
281 ax.loglog(T_range, eta, linewidth=2.5, color=color,
282 label=f'n_e = {n_e:.1e} m⁻³')
283
284 # Copper resistivity for reference
285 eta_copper = 1.7e-8 # Ω·m at room temperature
286 ax.axhline(eta_copper, color='brown', linestyle='--', linewidth=2,
287 label='Copper (room temp)')
288
289 ax.set_xlabel('Temperature [eV]', fontsize=13)
290 ax.set_ylabel('Resistivity η [Ω·m]', fontsize=13)
291 ax.set_title('Spitzer Resistivity vs Temperature',
292 fontsize=14, fontweight='bold')
293 ax.legend(loc='best', fontsize=10)
294 ax.grid(True, alpha=0.3, which='both')
295
296 plt.tight_layout()
297 plt.savefig('spitzer_resistivity.png', dpi=300, bbox_inches='tight')
298 plt.show()
299
300
301def plot_mean_free_path():
302 """Plot mean free path vs temperature."""
303
304 # Temperature range
305 T_range = np.logspace(0, 5, 200) # 1 eV to 100 keV
306
307 # Different densities
308 densities = [1e18, 1e19, 1e20, 1e21] # m^-3
309
310 fig, ax = plt.subplots(figsize=(10, 7))
311
312 colors = plt.cm.cool(np.linspace(0.1, 0.9, len(densities)))
313
314 for n_e, color in zip(densities, colors):
315 mfp = mean_free_path(n_e * np.ones_like(T_range), T_range)
316 ax.loglog(T_range, mfp, linewidth=2.5, color=color,
317 label=f'n_e = {n_e:.1e} m⁻³')
318
319 ax.set_xlabel('Temperature [eV]', fontsize=13)
320 ax.set_ylabel('Mean Free Path λ_mfp [m]', fontsize=13)
321 ax.set_title('Electron Mean Free Path vs Temperature',
322 fontsize=14, fontweight='bold')
323 ax.legend(loc='best', fontsize=10)
324 ax.grid(True, alpha=0.3, which='both')
325
326 # Add reference lines
327 ax.axhline(1.0, color='gray', linestyle=':', alpha=0.5, label='1 m')
328 ax.axhline(1e-3, color='gray', linestyle=':', alpha=0.5, label='1 mm')
329
330 plt.tight_layout()
331 plt.savefig('mean_free_path.png', dpi=300, bbox_inches='tight')
332 plt.show()
333
334
335def plot_timescale_comparison():
336 """Compare collision times with other plasma timescales."""
337
338 n_e = 1e20 # m^-3
339 B = 5.0 # Tesla
340 T_range = np.logspace(0, 5, 200) # 1 eV to 100 keV
341
342 # Collision time
343 nu_ei = electron_ion_collision_freq(n_e, T_range)
344 tau_coll = 1 / nu_ei
345
346 # Plasma period
347 omega_pe = np.sqrt(n_e * e**2 / (epsilon_0 * m_e))
348 tau_pe = 2 * np.pi / omega_pe
349
350 # Gyroperiod
351 omega_ce = e * B / m_e
352 tau_ce = 2 * np.pi / omega_ce
353
354 # Create figure
355 fig, ax = plt.subplots(figsize=(10, 8))
356
357 ax.loglog(T_range, tau_coll, 'r-', linewidth=3, label='τ_coll = 1/ν_ei')
358 ax.axhline(tau_pe, color='blue', linestyle='--', linewidth=2,
359 label=f'τ_pe = 2π/ω_pe = {tau_pe:.2e} s')
360 ax.axhline(tau_ce, color='green', linestyle='--', linewidth=2,
361 label=f'τ_ce = 2π/ω_ce = {tau_ce:.2e} s (B=5T)')
362
363 # Shaded regions
364 ax.fill_between(T_range, 1e-15, tau_pe, alpha=0.2, color='blue',
365 label='τ < τ_pe (plasma oscillation)')
366 ax.fill_between(T_range, tau_pe, tau_ce, alpha=0.2, color='green',
367 label='τ_pe < τ < τ_ce (gyration)')
368
369 ax.set_xlabel('Electron Temperature [eV]', fontsize=13)
370 ax.set_ylabel('Time [s]', fontsize=13)
371 ax.set_title(f'Plasma Timescales Comparison\n(n_e = {n_e:.1e} m⁻³, B = {B} T)',
372 fontsize=14, fontweight='bold')
373 ax.legend(loc='best', fontsize=10)
374 ax.grid(True, alpha=0.3, which='both')
375 ax.set_ylim([1e-15, 1e0])
376
377 # Add annotations
378 ax.text(10, 1e-6, 'Collisional', fontsize=12, color='red', fontweight='bold')
379 ax.text(1e4, 1e-6, 'Collisionless', fontsize=12, color='blue', fontweight='bold')
380
381 plt.tight_layout()
382 plt.savefig('timescale_comparison.png', dpi=300, bbox_inches='tight')
383 plt.show()
384
385
386def print_example_calculations():
387 """Print example calculations for different plasmas."""
388
389 print("\n" + "="*90)
390 print("COLLISION FREQUENCY CALCULATIONS FOR VARIOUS PLASMAS")
391 print("="*90 + "\n")
392
393 plasmas = [
394 ('Tokamak Core', 1e20, 10e3, 5.0),
395 ('Solar Corona', 1e14, 100, 1e-4),
396 ('Ionosphere', 1e11, 0.1, 50e-6),
397 ('Neon Sign', 1e16, 2, 0.0),
398 ]
399
400 for name, n_e, T_e, B in plasmas:
401 print(f"{name}:")
402 print(f" n_e = {n_e:.2e} m^-3, T_e = {T_e:.2e} eV, B = {B} T")
403
404 ln_Lambda = coulomb_logarithm(n_e, T_e)
405 nu_ei = electron_ion_collision_freq(n_e, T_e)
406 eta = spitzer_resistivity(n_e, T_e)
407 mfp = mean_free_path(n_e, T_e)
408
409 omega_pe = np.sqrt(n_e * e**2 / (epsilon_0 * m_e))
410 tau_pe = 1 / omega_pe
411
412 print(f" Coulomb logarithm: {ln_Lambda:.2f}")
413 print(f" Collision freq ν_ei: {nu_ei:.3e} Hz")
414 print(f" Collision time τ_coll: {1/nu_ei:.3e} s")
415 print(f" Plasma period τ_pe: {tau_pe:.3e} s")
416 print(f" Collisionality ν/ω_pe: {nu_ei * tau_pe:.3e}")
417 print(f" Spitzer resistivity: {eta:.3e} Ω·m")
418 print(f" Mean free path: {mfp:.3e} m")
419
420 if B > 0:
421 omega_ce = e * B / m_e
422 tau_ce = 1 / omega_ce
423 print(f" Gyroperiod τ_ce: {tau_ce:.3e} s")
424 print(f" ν_ei/ω_ce: {nu_ei / omega_ce:.3e}")
425
426 print()
427
428
429if __name__ == '__main__':
430 print("\n" + "="*90)
431 print("COLLISION FREQUENCIES AND TRANSPORT IN PLASMAS")
432 print("="*90)
433
434 # Print example calculations
435 print_example_calculations()
436
437 # Generate plots
438 print("Generating plots...")
439 print(" 1. Collision frequencies vs temperature...")
440 plot_collision_frequencies()
441
442 print(" 2. Coulomb logarithm...")
443 plot_coulomb_logarithm()
444
445 print(" 3. Spitzer resistivity...")
446 plot_spitzer_resistivity()
447
448 print(" 4. Mean free path...")
449 plot_mean_free_path()
450
451 print(" 5. Timescale comparison...")
452 plot_timescale_comparison()
453
454 print("\nDone! Generated 5 figures:")
455 print(" - collision_frequencies.png")
456 print(" - coulomb_logarithm.png")
457 print(" - spitzer_resistivity.png")
458 print(" - mean_free_path.png")
459 print(" - timescale_comparison.png")