1#!/usr/bin/env python3
2"""
3Ion Acoustic Wave Dispersion and Damping
4
5This script computes and visualizes ion acoustic wave properties including
6dispersion relation, damping rates, and the effect of electron-to-ion
7temperature ratio.
8
9Key Physics:
10- Sound speed: cs = sqrt(kTe/mi)
11- Dispersion: ω = k*cs / sqrt(1 + k²λD²)
12- Landau damping depends strongly on Te/Ti ratio
13- Weakly damped only when Te >> Ti
14
15Author: Plasma Physics Examples
16License: MIT
17"""
18
19import numpy as np
20import matplotlib.pyplot as plt
21from matplotlib.gridspec import GridSpec
22from scipy.special import erf
23
24# Physical constants
25EPS0 = 8.854187817e-12 # F/m
26QE = 1.602176634e-19 # C
27ME = 9.10938356e-31 # kg
28MP = 1.672621898e-27 # kg (proton mass)
29KB = 1.380649e-23 # J/K
30
31def compute_sound_speed(Te, Ti, mi):
32 """
33 Compute ion acoustic sound speed.
34
35 cs = sqrt((kTe + 3kTi) / mi) ≈ sqrt(kTe/mi) for Te >> Ti
36
37 Parameters:
38 -----------
39 Te : float
40 Electron temperature [eV]
41 Ti : float
42 Ion temperature [eV]
43 mi : float
44 Ion mass [kg]
45
46 Returns:
47 --------
48 cs : float
49 Sound speed [m/s]
50 """
51 Te_joule = Te * QE
52 Ti_joule = Ti * QE
53
54 # Full expression including ion pressure
55 cs = np.sqrt((Te_joule + 3 * Ti_joule) / mi)
56
57 return cs
58
59def ion_acoustic_dispersion(k, ne, Te, Ti, mi):
60 """
61 Compute ion acoustic wave dispersion relation.
62
63 ω = k*cs / sqrt(1 + k²λD²)
64
65 Parameters:
66 -----------
67 k : array
68 Wavenumber [rad/m]
69 ne : float
70 Electron density [m^-3]
71 Te : float
72 Electron temperature [eV]
73 Ti : float
74 Ion temperature [eV]
75 mi : float
76 Ion mass [kg]
77
78 Returns:
79 --------
80 omega : array
81 Angular frequency [rad/s]
82 """
83 # Compute sound speed
84 cs = compute_sound_speed(Te, Ti, mi)
85
86 # Debye length (using electron temperature)
87 Te_joule = Te * QE
88 lambda_D = np.sqrt(EPS0 * Te_joule / (ne * QE**2))
89
90 # Dispersion relation
91 omega = k * cs / np.sqrt(1 + k**2 * lambda_D**2)
92
93 return omega
94
95def ion_landau_damping(k, ne, Te, Ti, mi):
96 """
97 Compute ion Landau damping rate for ion acoustic waves.
98
99 This is an approximate formula valid for Te/Ti >> 1.
100
101 γ/ω ≈ -sqrt(π/8) * (me/mi)^(1/2) * (1 + Te/Ti)^(-3/2) * exp(-Te/(2Ti) - 3/2)
102
103 Parameters:
104 -----------
105 k : array
106 Wavenumber [rad/m]
107 ne : float
108 Electron density [m^-3]
109 Te : float
110 Electron temperature [eV]
111 Ti : float
112 Ion temperature [eV]
113 mi : float
114 Ion mass [kg]
115
116 Returns:
117 --------
118 gamma : array
119 Damping rate [rad/s]
120 """
121 # Compute frequency
122 omega = ion_acoustic_dispersion(k, ne, Te, Ti, mi)
123
124 # Temperature ratio
125 tau = Te / Ti
126
127 # Approximate damping rate (valid for tau >> 1)
128 # γ/ω ≈ -sqrt(π/8) * (1/tau)^(3/2) * exp(-1/(2*tau) - 3/2)
129 if tau > 1:
130 gamma_over_omega = -np.sqrt(np.pi / 8) * tau**(-3/2) * np.exp(-tau/2 - 3/2)
131 else:
132 # Strong damping when Te ~ Ti
133 gamma_over_omega = -1.0
134
135 gamma = gamma_over_omega * omega
136
137 return gamma
138
139def plot_ion_acoustic_waves():
140 """
141 Create comprehensive visualization of ion acoustic wave properties.
142 """
143 # Plasma parameters
144 ne = 1e18 # m^-3
145 mi = MP # Hydrogen ions
146
147 # Different temperature ratios to explore
148 tau_values = [1, 3, 10, 100] # Te/Ti
149 colors = ['red', 'orange', 'green', 'blue']
150
151 # Create wavenumber array
152 Te_ref = 10.0 # eV
153 lambda_D_ref = np.sqrt(EPS0 * Te_ref * QE / (ne * QE**2))
154 k = np.linspace(1e2, 10 / lambda_D_ref, 1000)
155
156 # Create figure
157 fig = plt.figure(figsize=(14, 10))
158 gs = GridSpec(3, 2, figure=fig, hspace=0.3, wspace=0.3)
159
160 # Store data for different tau values
161 omega_data = []
162 gamma_data = []
163 cs_data = []
164
165 print("=" * 70)
166 print("Ion Acoustic Wave Parameters")
167 print("=" * 70)
168 print(f"Electron density: {ne:.2e} m^-3")
169 print(f"Ion mass: {mi/MP:.1f} × proton mass")
170 print("-" * 70)
171
172 for tau, color in zip(tau_values, colors):
173 Te = 10.0 # eV
174 Ti = Te / tau
175
176 # Compute dispersion
177 omega = ion_acoustic_dispersion(k, ne, Te, Ti, mi)
178 gamma = ion_landau_damping(k, ne, Te, Ti, mi)
179 cs = compute_sound_speed(Te, Ti, mi)
180
181 omega_data.append(omega)
182 gamma_data.append(gamma)
183 cs_data.append(cs)
184
185 print(f"Te/Ti = {tau:6.1f}: Te = {Te:5.1f} eV, Ti = {Ti:5.2f} eV, "
186 f"cs = {cs/1e3:6.2f} km/s")
187
188 print("=" * 70)
189
190 # Plot 1: Dispersion relation ω(k)
191 ax1 = fig.add_subplot(gs[0, :])
192 for i, (tau, color) in enumerate(zip(tau_values, colors)):
193 k_lambda = k * lambda_D_ref
194 omega = omega_data[i]
195 cs = cs_data[i]
196
197 ax1.plot(k_lambda, omega / (k * cs),
198 color=color, linewidth=2, label=f'Te/Ti = {tau}')
199
200 # Plot long wavelength limit (should approach 1)
201 ax1.axhline(y=1.0, color='black', linestyle='--', linewidth=1,
202 label='Long wavelength limit')
203
204 ax1.set_xlabel(r'$k \lambda_D$', fontsize=12)
205 ax1.set_ylabel(r'$\omega / (k c_s)$', fontsize=12)
206 ax1.set_title('Ion Acoustic Wave Dispersion (Normalized)',
207 fontsize=14, fontweight='bold')
208 ax1.grid(True, alpha=0.3)
209 ax1.legend(fontsize=10, loc='upper right')
210 ax1.set_xlim([0, 10])
211 ax1.set_ylim([0, 1.1])
212
213 # Plot 2: Absolute frequency
214 ax2 = fig.add_subplot(gs[1, 0])
215 for i, (tau, color) in enumerate(zip(tau_values, colors)):
216 k_lambda = k * lambda_D_ref
217 omega = omega_data[i]
218
219 ax2.plot(k_lambda, omega / (2 * np.pi * 1e9),
220 color=color, linewidth=2, label=f'Te/Ti = {tau}')
221
222 ax2.set_xlabel(r'$k \lambda_D$', fontsize=12)
223 ax2.set_ylabel('Frequency (GHz)', fontsize=12)
224 ax2.set_title('Absolute Frequency', fontsize=12, fontweight='bold')
225 ax2.grid(True, alpha=0.3)
226 ax2.legend(fontsize=10)
227 ax2.set_xlim([0, 10])
228
229 # Plot 3: Phase velocity
230 ax3 = fig.add_subplot(gs[1, 1])
231 for i, (tau, color) in enumerate(zip(tau_values, colors)):
232 k_lambda = k * lambda_D_ref
233 omega = omega_data[i]
234 cs = cs_data[i]
235 vph = omega / k
236
237 ax3.plot(k_lambda, vph / cs,
238 color=color, linewidth=2, label=f'Te/Ti = {tau}')
239
240 ax3.axhline(y=1.0, color='black', linestyle='--', linewidth=1,
241 label='Sound speed')
242
243 ax3.set_xlabel(r'$k \lambda_D$', fontsize=12)
244 ax3.set_ylabel(r'$v_{ph} / c_s$', fontsize=12)
245 ax3.set_title('Phase Velocity', fontsize=12, fontweight='bold')
246 ax3.grid(True, alpha=0.3)
247 ax3.legend(fontsize=10)
248 ax3.set_xlim([0, 10])
249
250 # Plot 4: Damping rate (absolute)
251 ax4 = fig.add_subplot(gs[2, 0])
252 for i, (tau, color) in enumerate(zip(tau_values, colors)):
253 if tau >= 3: # Only plot for cases with weak damping
254 k_lambda = k * lambda_D_ref
255 gamma = gamma_data[i]
256
257 ax4.plot(k_lambda, -gamma / (2 * np.pi * 1e9),
258 color=color, linewidth=2, label=f'Te/Ti = {tau}')
259
260 ax4.set_xlabel(r'$k \lambda_D$', fontsize=12)
261 ax4.set_ylabel('Damping Rate (GHz)', fontsize=12)
262 ax4.set_title('Ion Landau Damping Rate (Te/Ti ≥ 3)',
263 fontsize=12, fontweight='bold')
264 ax4.grid(True, alpha=0.3)
265 ax4.legend(fontsize=10)
266 ax4.set_xlim([0, 10])
267 ax4.set_yscale('log')
268
269 # Plot 5: Damping rate normalized
270 ax5 = fig.add_subplot(gs[2, 1])
271 for i, (tau, color) in enumerate(zip(tau_values, colors)):
272 if tau >= 3: # Only plot for cases with weak damping
273 k_lambda = k * lambda_D_ref
274 omega = omega_data[i]
275 gamma = gamma_data[i]
276
277 ax5.plot(k_lambda, -gamma / omega,
278 color=color, linewidth=2, label=f'Te/Ti = {tau}')
279
280 ax5.set_xlabel(r'$k \lambda_D$', fontsize=12)
281 ax5.set_ylabel(r'$|\gamma| / \omega$', fontsize=12)
282 ax5.set_title('Normalized Damping Rate (Te/Ti ≥ 3)',
283 fontsize=12, fontweight='bold')
284 ax5.grid(True, alpha=0.3)
285 ax5.legend(fontsize=10)
286 ax5.set_xlim([0, 10])
287 ax5.set_yscale('log')
288
289 plt.suptitle('Ion Acoustic Wave: Effect of Temperature Ratio Te/Ti',
290 fontsize=16, fontweight='bold', y=0.995)
291
292 plt.savefig('ion_acoustic_wave.png', dpi=150, bbox_inches='tight')
293 print("\nPlot saved as 'ion_acoustic_wave.png'")
294
295 plt.show()
296
297if __name__ == "__main__":
298 plot_ion_acoustic_waves()