1#!/usr/bin/env python3
2"""
3Debye Shielding Visualization
4
5This script demonstrates the shielding of a test charge in plasma by solving
6the Poisson equation with Boltzmann electrons. Shows both linearized analytical
7solution and numerical nonlinear solution.
8
9Author: Plasma Physics Examples
10License: MIT
11"""
12
13import numpy as np
14import matplotlib.pyplot as plt
15from scipy.constants import e, epsilon_0, k
16from scipy.integrate import odeint
17
18
19def debye_length(n0, T_eV):
20 """
21 Calculate Debye length.
22
23 Parameters:
24 -----------
25 n0 : float
26 Background density [m^-3]
27 T_eV : float
28 Temperature [eV]
29
30 Returns:
31 --------
32 float : Debye length [m]
33 """
34 T_J = T_eV * e
35 return np.sqrt(epsilon_0 * T_J / (n0 * e**2))
36
37
38def bare_coulomb_potential(r, q):
39 """
40 Bare Coulomb potential.
41
42 Parameters:
43 -----------
44 r : array
45 Radial distance [m]
46 q : float
47 Test charge [C]
48
49 Returns:
50 --------
51 array : Potential [V]
52 """
53 # Avoid singularity at r=0
54 r_safe = np.where(r > 0, r, 1e-10)
55 return q / (4 * np.pi * epsilon_0 * r_safe)
56
57
58def debye_shielded_potential(r, q, lambda_D):
59 """
60 Debye-shielded potential (analytical, linearized).
61
62 Parameters:
63 -----------
64 r : array
65 Radial distance [m]
66 q : float
67 Test charge [C]
68 lambda_D : float
69 Debye length [m]
70
71 Returns:
72 --------
73 array : Potential [V]
74 """
75 r_safe = np.where(r > 0, r, 1e-10)
76 return (q / (4 * np.pi * epsilon_0 * r_safe)) * np.exp(-r_safe / lambda_D)
77
78
79def poisson_1d_spherical(y, r, n0, T_eV, q):
80 """
81 1D spherical Poisson equation with Boltzmann electrons.
82
83 d²φ/dr² + (2/r)dφ/dr = (e*n0/ε0)(exp(eφ/kT) - 1)
84
85 Convert to system of ODEs:
86 y[0] = φ
87 y[1] = dφ/dr
88
89 Parameters:
90 -----------
91 y : array
92 [φ, dφ/dr]
93 r : float
94 Radial coordinate [m]
95 n0 : float
96 Background density [m^-3]
97 T_eV : float
98 Temperature [eV]
99 q : float
100 Test charge [C]
101
102 Returns:
103 --------
104 array : [dφ/dr, d²φ/dr²]
105 """
106 phi, dphi_dr = y
107 T_J = T_eV * e
108
109 # Avoid singularity at r=0
110 if r < 1e-12:
111 r = 1e-12
112
113 # Poisson equation
114 d2phi_dr2 = (e * n0 / epsilon_0) * (np.exp(e * phi / T_J) - 1) - (2 / r) * dphi_dr
115
116 return [dphi_dr, d2phi_dr2]
117
118
119def solve_poisson_numerical(r_array, n0, T_eV, q):
120 """
121 Solve Poisson equation numerically with boundary conditions.
122
123 Parameters:
124 -----------
125 r_array : array
126 Radial grid [m]
127 n0 : float
128 Background density [m^-3]
129 T_eV : float
130 Temperature [eV]
131 q : float
132 Test charge [C]
133
134 Returns:
135 --------
136 array : Potential φ(r) [V]
137 """
138 # Boundary condition at small r (approximately bare Coulomb)
139 r0 = r_array[0]
140 phi0 = q / (4 * np.pi * epsilon_0 * r0)
141 dphi0 = -q / (4 * np.pi * epsilon_0 * r0**2)
142
143 # Initial condition
144 y0 = [phi0, dphi0]
145
146 # Solve ODE
147 solution = odeint(poisson_1d_spherical, y0, r_array, args=(n0, T_eV, q))
148
149 return solution[:, 0]
150
151
152def plot_shielding_comparison():
153 """Plot comparison of bare Coulomb, linearized Debye, and numerical solution."""
154
155 # Parameters
156 n0 = 1e16 # m^-3
157 T_eV = 2.0 # eV
158 q = e # Test charge (one elementary charge)
159
160 lambda_D = debye_length(n0, T_eV)
161
162 # Radial grid
163 r = np.logspace(-6, -2, 500) # 1 μm to 1 cm
164
165 # Compute potentials
166 phi_bare = bare_coulomb_potential(r, q)
167 phi_debye = debye_shielded_potential(r, q, lambda_D)
168 phi_numerical = solve_poisson_numerical(r, n0, T_eV, q)
169
170 # Create figure
171 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
172
173 # Plot 1: Linear scale
174 ax1.plot(r * 1e3, phi_bare, 'k--', linewidth=2, label='Bare Coulomb', alpha=0.7)
175 ax1.plot(r * 1e3, phi_debye, 'b-', linewidth=2, label='Debye (linearized)')
176 ax1.plot(r * 1e3, phi_numerical, 'r:', linewidth=2.5, label='Numerical (nonlinear)')
177 ax1.axvline(lambda_D * 1e3, color='gray', linestyle='--', alpha=0.5,
178 label=f'λ_D = {lambda_D*1e3:.3f} mm')
179
180 ax1.set_xlabel('Distance r [mm]', fontsize=12)
181 ax1.set_ylabel('Potential φ(r) [V]', fontsize=12)
182 ax1.set_title('Debye Shielding of Test Charge', fontsize=13, fontweight='bold')
183 ax1.legend(loc='best', fontsize=10)
184 ax1.grid(True, alpha=0.3)
185 ax1.set_xlim([0, 5])
186 ax1.set_ylim([0, 20])
187
188 # Plot 2: Log-log scale
189 ax2.loglog(r * 1e3, phi_bare, 'k--', linewidth=2, label='Bare Coulomb', alpha=0.7)
190 ax2.loglog(r * 1e3, phi_debye, 'b-', linewidth=2, label='Debye (linearized)')
191 ax2.loglog(r * 1e3, phi_numerical, 'r:', linewidth=2.5, label='Numerical (nonlinear)')
192 ax2.axvline(lambda_D * 1e3, color='gray', linestyle='--', alpha=0.5,
193 label=f'λ_D = {lambda_D*1e3:.3f} mm')
194
195 ax2.set_xlabel('Distance r [mm]', fontsize=12)
196 ax2.set_ylabel('Potential φ(r) [V]', fontsize=12)
197 ax2.set_title('Debye Shielding (Log Scale)', fontsize=13, fontweight='bold')
198 ax2.legend(loc='best', fontsize=10)
199 ax2.grid(True, alpha=0.3, which='both')
200
201 plt.tight_layout()
202 plt.savefig('debye_shielding_comparison.png', dpi=300, bbox_inches='tight')
203 plt.show()
204
205
206def plot_temperature_dependence():
207 """Show how shielding changes with temperature."""
208
209 # Parameters
210 n0 = 1e16 # m^-3
211 q = e
212 temperatures = [0.5, 1.0, 2.0, 5.0] # eV
213
214 # Radial grid
215 r = np.linspace(0.001e-3, 5e-3, 200) # mm scale
216
217 # Create figure
218 fig, ax = plt.subplots(figsize=(10, 7))
219
220 colors = plt.cm.viridis(np.linspace(0.2, 0.9, len(temperatures)))
221
222 for T_eV, color in zip(temperatures, colors):
223 lambda_D = debye_length(n0, T_eV)
224 phi_debye = debye_shielded_potential(r, q, lambda_D)
225
226 ax.plot(r * 1e3, phi_debye, linewidth=2.5, color=color,
227 label=f'T = {T_eV} eV, λ_D = {lambda_D*1e3:.3f} mm')
228 ax.axvline(lambda_D * 1e3, color=color, linestyle=':', alpha=0.5)
229
230 # Bare Coulomb for reference
231 phi_bare = bare_coulomb_potential(r, q)
232 ax.plot(r * 1e3, phi_bare, 'k--', linewidth=2, label='Bare Coulomb', alpha=0.5)
233
234 ax.set_xlabel('Distance r [mm]', fontsize=12)
235 ax.set_ylabel('Potential φ(r) [V]', fontsize=12)
236 ax.set_title(f'Temperature Dependence of Debye Shielding\n(n₀ = {n0:.1e} m⁻³)',
237 fontsize=13, fontweight='bold')
238 ax.legend(loc='best', fontsize=10)
239 ax.grid(True, alpha=0.3)
240 ax.set_xlim([0, 5])
241 ax.set_ylim([0, 30])
242
243 plt.tight_layout()
244 plt.savefig('debye_temperature_dependence.png', dpi=300, bbox_inches='tight')
245 plt.show()
246
247
248def plot_density_response():
249 """Visualize electron density response around test charge."""
250
251 # Parameters
252 n0 = 1e16 # m^-3
253 T_eV = 2.0 # eV
254 q = e
255
256 lambda_D = debye_length(n0, T_eV)
257 T_J = T_eV * e
258
259 # Radial grid
260 r = np.linspace(0.001e-3, 10e-3, 500) # mm scale
261
262 # Compute potential
263 phi_debye = debye_shielded_potential(r, q, lambda_D)
264
265 # Electron density from Boltzmann relation
266 n_e = n0 * np.exp(e * phi_debye / T_J)
267
268 # Create figure
269 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 10))
270
271 # Plot 1: Potential
272 ax1.plot(r * 1e3, phi_debye, 'b-', linewidth=2.5)
273 ax1.axvline(lambda_D * 1e3, color='red', linestyle='--', linewidth=2,
274 label=f'λ_D = {lambda_D*1e3:.3f} mm')
275 ax1.axhline(0, color='k', linestyle='-', linewidth=0.8)
276
277 ax1.set_xlabel('Distance r [mm]', fontsize=12)
278 ax1.set_ylabel('Potential φ(r) [V]', fontsize=12)
279 ax1.set_title('Electrostatic Potential', fontsize=13, fontweight='bold')
280 ax1.legend(loc='best', fontsize=11)
281 ax1.grid(True, alpha=0.3)
282
283 # Plot 2: Density perturbation
284 delta_n = n_e - n0
285 ax2.plot(r * 1e3, delta_n / n0 * 100, 'r-', linewidth=2.5)
286 ax2.axvline(lambda_D * 1e3, color='red', linestyle='--', linewidth=2,
287 label=f'λ_D = {lambda_D*1e3:.3f} mm')
288 ax2.axhline(0, color='k', linestyle='-', linewidth=0.8)
289
290 ax2.set_xlabel('Distance r [mm]', fontsize=12)
291 ax2.set_ylabel('Density Perturbation Δn/n₀ [%]', fontsize=12)
292 ax2.set_title('Electron Density Response', fontsize=13, fontweight='bold')
293 ax2.legend(loc='best', fontsize=11)
294 ax2.grid(True, alpha=0.3)
295
296 plt.tight_layout()
297 plt.savefig('debye_density_response.png', dpi=300, bbox_inches='tight')
298 plt.show()
299
300
301def plot_2d_potential():
302 """2D visualization of shielded potential."""
303
304 # Parameters
305 n0 = 1e16 # m^-3
306 T_eV = 2.0 # eV
307 q = e
308
309 lambda_D = debye_length(n0, T_eV)
310
311 # Create 2D grid
312 x = np.linspace(-5e-3, 5e-3, 200)
313 y = np.linspace(-5e-3, 5e-3, 200)
314 X, Y = np.meshgrid(x, y)
315 R = np.sqrt(X**2 + Y**2)
316
317 # Compute potentials
318 phi_bare = bare_coulomb_potential(R, q)
319 phi_debye = debye_shielded_potential(R, q, lambda_D)
320
321 # Create figure
322 fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 5))
323
324 # Plot 1: Bare Coulomb
325 levels = np.linspace(0, 20, 30)
326 cs1 = ax1.contourf(X * 1e3, Y * 1e3, phi_bare, levels=levels, cmap='hot')
327 ax1.contour(X * 1e3, Y * 1e3, phi_bare, levels=[1, 5, 10, 15],
328 colors='white', linewidths=1, alpha=0.5)
329 ax1.set_xlabel('x [mm]', fontsize=12)
330 ax1.set_ylabel('y [mm]', fontsize=12)
331 ax1.set_title('Bare Coulomb Potential', fontsize=13, fontweight='bold')
332 ax1.set_aspect('equal')
333 plt.colorbar(cs1, ax=ax1, label='φ [V]')
334
335 # Plot 2: Debye shielded
336 cs2 = ax2.contourf(X * 1e3, Y * 1e3, phi_debye, levels=levels, cmap='hot')
337 ax2.contour(X * 1e3, Y * 1e3, phi_debye, levels=[1, 5, 10, 15],
338 colors='white', linewidths=1, alpha=0.5)
339 circle = plt.Circle((0, 0), lambda_D * 1e3, color='cyan', fill=False,
340 linewidth=2, linestyle='--', label='λ_D')
341 ax2.add_patch(circle)
342 ax2.set_xlabel('x [mm]', fontsize=12)
343 ax2.set_ylabel('y [mm]', fontsize=12)
344 ax2.set_title('Debye-Shielded Potential', fontsize=13, fontweight='bold')
345 ax2.set_aspect('equal')
346 ax2.legend(loc='upper right', fontsize=10)
347 plt.colorbar(cs2, ax=ax2, label='φ [V]')
348
349 # Plot 3: Difference
350 diff = phi_bare - phi_debye
351 cs3 = ax3.contourf(X * 1e3, Y * 1e3, diff, levels=20, cmap='coolwarm')
352 ax3.contour(X * 1e3, Y * 1e3, diff, levels=10, colors='black',
353 linewidths=0.5, alpha=0.3)
354 circle = plt.Circle((0, 0), lambda_D * 1e3, color='green', fill=False,
355 linewidth=2, linestyle='--', label='λ_D')
356 ax3.add_patch(circle)
357 ax3.set_xlabel('x [mm]', fontsize=12)
358 ax3.set_ylabel('y [mm]', fontsize=12)
359 ax3.set_title('Shielding Effect (Difference)', fontsize=13, fontweight='bold')
360 ax3.set_aspect('equal')
361 ax3.legend(loc='upper right', fontsize=10)
362 plt.colorbar(cs3, ax=ax3, label='Δφ [V]')
363
364 plt.tight_layout()
365 plt.savefig('debye_2d_potential.png', dpi=300, bbox_inches='tight')
366 plt.show()
367
368
369if __name__ == '__main__':
370 print("\n" + "="*80)
371 print("DEBYE SHIELDING VISUALIZATION")
372 print("="*80 + "\n")
373
374 # Example parameters
375 n0 = 1e16 # m^-3
376 T_eV = 2.0 # eV
377 q = e
378
379 lambda_D = debye_length(n0, T_eV)
380
381 print(f"Plasma parameters:")
382 print(f" Density: {n0:.2e} m^-3")
383 print(f" Temperature: {T_eV} eV")
384 print(f" Test charge: {q/e:.1f} e")
385 print(f"\nDebye length: {lambda_D:.4e} m = {lambda_D*1e3:.4f} mm\n")
386
387 print("Generating plots...")
388 print(" 1. Shielding comparison (linear vs numerical)...")
389 plot_shielding_comparison()
390
391 print(" 2. Temperature dependence...")
392 plot_temperature_dependence()
393
394 print(" 3. Density response...")
395 plot_density_response()
396
397 print(" 4. 2D potential visualization...")
398 plot_2d_potential()
399
400 print("\nDone! Generated 4 figures:")
401 print(" - debye_shielding_comparison.png")
402 print(" - debye_temperature_dependence.png")
403 print(" - debye_density_response.png")
404 print(" - debye_2d_potential.png")