1#!/usr/bin/env python3
2"""
3E×B Drift Simulation
4
5This script demonstrates the E×B drift of charged particles in crossed
6electric and magnetic fields. Shows that the drift is charge-independent
7and verifies the drift velocity formula.
8
9Author: Plasma Physics Examples
10License: MIT
11"""
12
13import numpy as np
14import matplotlib.pyplot as plt
15from scipy.constants import e, m_e, m_p
16
17
18def boris_push(x, v, q, m, E, B, dt):
19 """
20 Boris algorithm for particle pushing in electromagnetic fields.
21
22 Parameters:
23 -----------
24 x : array (3,)
25 Position [m]
26 v : array (3,)
27 Velocity [m/s]
28 q : float
29 Charge [C]
30 m : float
31 Mass [kg]
32 E : array (3,)
33 Electric field [V/m]
34 B : array (3,)
35 Magnetic field [T]
36 dt : float
37 Timestep [s]
38
39 Returns:
40 --------
41 x_new, v_new : Updated position and velocity
42 """
43 # Half acceleration from E field
44 v_minus = v + (q * E / m) * (dt / 2)
45
46 # Rotation from B field
47 t = (q * dt / (2 * m)) * B
48 t_mag2 = np.dot(t, t)
49 s = 2 * t / (1 + t_mag2)
50
51 v_prime = v_minus + np.cross(v_minus, t)
52 v_plus = v_minus + np.cross(v_prime, s)
53
54 # Half acceleration from E field
55 v_new = v_plus + (q * E / m) * (dt / 2)
56
57 # Update position
58 x_new = x + v_new * dt
59
60 return x_new, v_new
61
62
63def simulate_particle(q, m, E, B, v0, x0, dt, n_steps):
64 """
65 Simulate particle motion in E and B fields.
66
67 Parameters:
68 -----------
69 q : float
70 Charge [C]
71 m : float
72 Mass [kg]
73 E : array (3,)
74 Electric field [V/m]
75 B : array (3,)
76 Magnetic field [T]
77 v0 : array (3,)
78 Initial velocity [m/s]
79 x0 : array (3,)
80 Initial position [m]
81 dt : float
82 Timestep [s]
83 n_steps : int
84 Number of steps
85
86 Returns:
87 --------
88 trajectory : dict with 'x', 'y', 'z', 'vx', 'vy', 'vz', 't'
89 """
90 # Initialize arrays
91 x = np.zeros((n_steps, 3))
92 v = np.zeros((n_steps, 3))
93 t = np.zeros(n_steps)
94
95 x[0] = x0
96 v[0] = v0
97
98 # Time integration
99 for i in range(n_steps - 1):
100 x[i+1], v[i+1] = boris_push(x[i], v[i], q, m, E, B, dt)
101 t[i+1] = t[i] + dt
102
103 return {
104 'x': x[:, 0], 'y': x[:, 1], 'z': x[:, 2],
105 'vx': v[:, 0], 'vy': v[:, 1], 'vz': v[:, 2],
106 't': t
107 }
108
109
110def plot_exb_drift():
111 """Demonstrate E×B drift for electron and ion."""
112
113 # Fields configuration
114 E0 = 1000 # V/m
115 B0 = 0.1 # T
116 E = np.array([E0, 0, 0]) # E in x direction
117 B = np.array([0, 0, B0]) # B in z direction
118
119 # Theoretical E×B drift velocity
120 v_ExB = np.cross(E, B) / B0**2
121 v_ExB_mag = np.linalg.norm(v_ExB)
122
123 print(f"E = {E0} V/m (x direction)")
124 print(f"B = {B0} T (z direction)")
125 print(f"E×B drift velocity: {v_ExB} m/s")
126 print(f"Magnitude: {v_ExB_mag} m/s = {v_ExB_mag/1e3:.1f} km/s\n")
127
128 # Particle parameters
129 q_e = -e
130 q_p = e
131 m_e_kg = m_e
132 m_p_kg = m_p
133
134 # Initial conditions (small perpendicular velocity)
135 v0 = np.array([0, 1e4, 0]) # Small v_y
136 x0_e = np.array([0, 0, 0])
137 x0_p = np.array([0, 1, 0]) # Offset for visibility
138
139 # Time parameters
140 omega_ce = abs(q_e) * B0 / m_e_kg
141 T_ce = 2 * np.pi / omega_ce
142 dt = T_ce / 100
143 n_steps = 2000
144
145 # Simulate
146 traj_e = simulate_particle(q_e, m_e_kg, E, B, v0, x0_e, dt, n_steps)
147 traj_p = simulate_particle(q_p, m_p_kg, E, B, v0, x0_p, dt, n_steps)
148
149 # Calculate drift velocities from simulation
150 # Use linear fit to position vs time
151 t_start = int(0.2 * n_steps) # Skip initial transient
152 t_end = n_steps
153
154 drift_e_y = np.polyfit(traj_e['t'][t_start:t_end],
155 traj_e['y'][t_start:t_end], 1)[0]
156 drift_p_y = np.polyfit(traj_p['t'][t_start:t_end],
157 traj_p['y'][t_start:t_end], 1)[0]
158
159 print(f"Simulated drift velocities:")
160 print(f" Electron: {drift_e_y:.2f} m/s (y direction)")
161 print(f" Proton: {drift_p_y:.2f} m/s (y direction)")
162 print(f" Theory: {v_ExB[1]:.2f} m/s (y direction)")
163 print(f" Error (electron): {abs(drift_e_y - v_ExB[1])/v_ExB[1] * 100:.2f}%")
164 print(f" Error (proton): {abs(drift_p_y - v_ExB[1])/v_ExB[1] * 100:.2f}%\n")
165
166 # Create figure
167 fig, axes = plt.subplots(2, 2, figsize=(14, 12))
168
169 # Plot 1: Trajectories in x-y plane
170 ax1 = axes[0, 0]
171 ax1.plot(traj_e['x'], traj_e['y'], 'b-', linewidth=1.5, label='Electron', alpha=0.7)
172 ax1.plot(traj_p['x'], traj_p['y'], 'r-', linewidth=1.5, label='Proton', alpha=0.7)
173 ax1.plot(traj_e['x'][0], traj_e['y'][0], 'go', markersize=10, label='Start (e⁻)')
174 ax1.plot(traj_p['x'][0], traj_p['y'][0], 'mo', markersize=10, label='Start (p⁺)')
175
176 # Draw drift velocity arrow
177 y_mid = (traj_e['y'][0] + traj_p['y'][0]) / 2
178 ax1.arrow(0, y_mid, 0, 0.3, head_width=0.05, head_length=0.05,
179 fc='green', ec='green', linewidth=2)
180 ax1.text(0.1, y_mid + 0.15, 'v_ExB', fontsize=12, fontweight='bold', color='green')
181
182 # Draw field vectors
183 ax1.arrow(-0.3, 0.5, 0.2, 0, head_width=0.05, head_length=0.03,
184 fc='red', ec='red', linewidth=2)
185 ax1.text(-0.15, 0.6, 'E', fontsize=12, fontweight='bold', color='red')
186
187 ax1.set_xlabel('x [m]', fontsize=12)
188 ax1.set_ylabel('y [m]', fontsize=12)
189 ax1.set_title('E×B Drift Trajectories', fontsize=13, fontweight='bold')
190 ax1.legend(loc='best', fontsize=10)
191 ax1.grid(True, alpha=0.3)
192 ax1.set_aspect('equal')
193
194 # Plot 2: y position vs time
195 ax2 = axes[0, 1]
196 ax2.plot(traj_e['t'] * 1e6, traj_e['y'], 'b-', linewidth=2, label='Electron')
197 ax2.plot(traj_p['t'] * 1e6, traj_p['y'], 'r-', linewidth=2, label='Proton')
198
199 # Plot linear fit lines
200 t_fit = traj_e['t'][t_start:t_end]
201 ax2.plot(t_fit * 1e6, drift_e_y * t_fit + traj_e['y'][t_start],
202 'b--', linewidth=2, alpha=0.5, label=f'Fit: {drift_e_y:.1f} m/s')
203 ax2.plot(t_fit * 1e6, drift_p_y * t_fit + traj_p['y'][t_start],
204 'r--', linewidth=2, alpha=0.5, label=f'Fit: {drift_p_y:.1f} m/s')
205
206 ax2.set_xlabel('Time [μs]', fontsize=12)
207 ax2.set_ylabel('y Position [m]', fontsize=12)
208 ax2.set_title('Drift in y Direction (Both Species)', fontsize=13, fontweight='bold')
209 ax2.legend(loc='best', fontsize=10)
210 ax2.grid(True, alpha=0.3)
211
212 # Plot 3: Velocity components (electron)
213 ax3 = axes[1, 0]
214 ax3.plot(traj_e['t'] * 1e6, traj_e['vx'] / 1e3, 'b-', linewidth=2, label='v_x')
215 ax3.plot(traj_e['t'] * 1e6, traj_e['vy'] / 1e3, 'r-', linewidth=2, label='v_y')
216 ax3.axhline(v_ExB[1] / 1e3, color='green', linestyle='--', linewidth=2,
217 label=f'v_ExB = {v_ExB[1]/1e3:.1f} km/s')
218
219 ax3.set_xlabel('Time [μs]', fontsize=12)
220 ax3.set_ylabel('Velocity [km/s]', fontsize=12)
221 ax3.set_title('Electron Velocity Components', fontsize=13, fontweight='bold')
222 ax3.legend(loc='best', fontsize=10)
223 ax3.grid(True, alpha=0.3)
224
225 # Plot 4: Phase space (v_x vs v_y)
226 ax4 = axes[1, 1]
227 ax4.plot(traj_e['vx'] / 1e3, traj_e['vy'] / 1e3, 'b-',
228 linewidth=1, alpha=0.5, label='Electron')
229 ax4.plot(traj_p['vx'] / 1e3, traj_p['vy'] / 1e3, 'r-',
230 linewidth=1, alpha=0.5, label='Proton')
231 ax4.plot(0, v_ExB[1] / 1e3, 'g*', markersize=20,
232 label='E×B drift velocity')
233
234 ax4.set_xlabel('v_x [km/s]', fontsize=12)
235 ax4.set_ylabel('v_y [km/s]', fontsize=12)
236 ax4.set_title('Velocity Phase Space', fontsize=13, fontweight='bold')
237 ax4.legend(loc='best', fontsize=10)
238 ax4.grid(True, alpha=0.3)
239
240 plt.tight_layout()
241 plt.savefig('exb_drift.png', dpi=300, bbox_inches='tight')
242 plt.show()
243
244
245def plot_parallel_acceleration():
246 """Demonstrate acceleration along B when E is parallel to B."""
247
248 # Fields configuration: E parallel to B
249 E0 = 1000 # V/m
250 B0 = 0.1 # T
251 E = np.array([0, 0, E0]) # E in z direction
252 B = np.array([0, 0, B0]) # B in z direction
253
254 print(f"\nParallel E field case:")
255 print(f"E = {E0} V/m (z direction, parallel to B)")
256 print(f"B = {B0} T (z direction)")
257
258 # Particle parameters
259 q_e = -e
260 q_p = e
261 m_e_kg = m_e
262 m_p_kg = m_p
263
264 # Initial conditions
265 v0 = np.array([1e5, 0, 0]) # Initial perpendicular velocity
266 x0_e = np.array([0, 0, 0])
267 x0_p = np.array([0, 0, 0])
268
269 # Time parameters
270 omega_ce = abs(q_e) * B0 / m_e_kg
271 T_ce = 2 * np.pi / omega_ce
272 dt = T_ce / 100
273 n_steps = 1000
274
275 # Simulate
276 traj_e = simulate_particle(q_e, m_e_kg, E, B, v0, x0_e, dt, n_steps)
277 traj_p = simulate_particle(q_p, m_p_kg, E, B, v0, x0_p, dt, n_steps)
278
279 # Expected acceleration
280 a_e = q_e * E0 / m_e_kg
281 a_p = q_p * E0 / m_p_kg
282
283 # Create figure
284 fig, axes = plt.subplots(2, 2, figsize=(14, 12))
285
286 # Plot 1: 3D helical trajectory with acceleration
287 from mpl_toolkits.mplot3d import Axes3D
288 ax1 = plt.subplot(2, 2, 1, projection='3d')
289
290 ax1.plot(traj_e['x'] * 1e3, traj_e['y'] * 1e3, traj_e['z'] * 1e3,
291 'b-', linewidth=2, label='Electron')
292 ax1.plot([traj_e['x'][0] * 1e3], [traj_e['y'][0] * 1e3], [traj_e['z'][0] * 1e3],
293 'go', markersize=10, label='Start')
294
295 ax1.set_xlabel('x [mm]', fontsize=11)
296 ax1.set_ylabel('y [mm]', fontsize=11)
297 ax1.set_zlabel('z [mm]', fontsize=11)
298 ax1.set_title('Electron Trajectory (E ∥ B)', fontsize=13, fontweight='bold')
299 ax1.legend(loc='best', fontsize=9)
300
301 # Plot 2: Parallel velocity vs time
302 ax2 = axes[0, 1]
303 ax2.plot(traj_e['t'] * 1e6, traj_e['vz'] / 1e3, 'b-', linewidth=2, label='Electron')
304 ax2.plot(traj_p['t'] * 1e6, traj_p['vz'] / 1e3, 'r-', linewidth=2, label='Proton')
305
306 # Expected from kinematics: v = v0 + at
307 ax2.plot(traj_e['t'] * 1e6, (v0[2] + a_e * traj_e['t']) / 1e3,
308 'b--', linewidth=2, alpha=0.5, label='Theory (e⁻)')
309 ax2.plot(traj_p['t'] * 1e6, (v0[2] + a_p * traj_p['t']) / 1e3,
310 'r--', linewidth=2, alpha=0.5, label='Theory (p⁺)')
311
312 ax2.set_xlabel('Time [μs]', fontsize=12)
313 ax2.set_ylabel('v_z (parallel) [km/s]', fontsize=12)
314 ax2.set_title('Parallel Acceleration', fontsize=13, fontweight='bold')
315 ax2.legend(loc='best', fontsize=10)
316 ax2.grid(True, alpha=0.3)
317
318 # Plot 3: Perpendicular velocity (should remain constant in magnitude)
319 ax3 = axes[1, 0]
320 v_perp_e = np.sqrt(traj_e['vx']**2 + traj_e['vy']**2)
321 v_perp_p = np.sqrt(traj_p['vx']**2 + traj_p['vy']**2)
322
323 ax3.plot(traj_e['t'] * 1e6, v_perp_e / 1e3, 'b-', linewidth=2, label='Electron')
324 ax3.plot(traj_p['t'] * 1e6, v_perp_p / 1e3, 'r-', linewidth=2, label='Proton')
325 ax3.axhline(np.linalg.norm(v0) / 1e3, color='green', linestyle='--',
326 linewidth=2, label='Initial v_⊥')
327
328 ax3.set_xlabel('Time [μs]', fontsize=12)
329 ax3.set_ylabel('|v_⊥| [km/s]', fontsize=12)
330 ax3.set_title('Perpendicular Velocity (Conserved)', fontsize=13, fontweight='bold')
331 ax3.legend(loc='best', fontsize=10)
332 ax3.grid(True, alpha=0.3)
333
334 # Plot 4: Total kinetic energy
335 ax4 = axes[1, 1]
336 KE_e = 0.5 * m_e_kg * (traj_e['vx']**2 + traj_e['vy']**2 + traj_e['vz']**2)
337 KE_p = 0.5 * m_p_kg * (traj_p['vx']**2 + traj_p['vy']**2 + traj_p['vz']**2)
338
339 ax4.plot(traj_e['t'] * 1e6, KE_e / e, 'b-', linewidth=2, label='Electron')
340 ax4.plot(traj_p['t'] * 1e6, KE_p / e, 'r-', linewidth=2, label='Proton')
341
342 ax4.set_xlabel('Time [μs]', fontsize=12)
343 ax4.set_ylabel('Kinetic Energy [eV]', fontsize=12)
344 ax4.set_title('Energy Gain from E Field', fontsize=13, fontweight='bold')
345 ax4.legend(loc='best', fontsize=10)
346 ax4.grid(True, alpha=0.3)
347
348 plt.tight_layout()
349 plt.savefig('parallel_acceleration.png', dpi=300, bbox_inches='tight')
350 plt.show()
351
352
353if __name__ == '__main__':
354 print("\n" + "="*80)
355 print("E×B DRIFT SIMULATION")
356 print("="*80 + "\n")
357
358 print("Part 1: E×B Drift (E ⊥ B)")
359 print("-" * 80)
360 plot_exb_drift()
361
362 print("\nPart 2: Parallel Acceleration (E ∥ B)")
363 print("-" * 80)
364 plot_parallel_acceleration()
365
366 print("\nDone! Generated 2 figures:")
367 print(" - exb_drift.png")
368 print(" - parallel_acceleration.png")