1#!/usr/bin/env python3
2"""
3Alfvén Wave Collision and Energy Cascade
4
5This script simulates the collision of counter-propagating Alfvén wave packets
6and demonstrates how wave-wave interactions lead to energy cascade in MHD
7turbulence. The simulation illustrates:
8- Critical balance: perpendicular cascade time ~ Alfvén wave period
9- Energy transfer from large to small scales
10- Formation of smaller-scale structures via nonlinear interactions
11
12Key results:
13- Counter-propagating waves generate small-scale perturbations
14- Energy cascades to higher wavenumbers
15- Critical balance τ_nl ~ τ_A determines cascade dynamics
16
17Author: Claude
18Date: 2026-02-14
19"""
20
21import numpy as np
22import matplotlib.pyplot as plt
23from matplotlib.animation import FuncAnimation
24from scipy.fft import fft, ifft, fftfreq
25
26
27def alfven_wave_packet(x, t, k0, width, phase, amplitude, v_alfven, direction=1):
28 """
29 Generate an Alfvén wave packet.
30
31 Parameters
32 ----------
33 x : ndarray
34 Spatial grid
35 t : float
36 Time
37 k0 : float
38 Central wavenumber
39 width : float
40 Packet width
41 phase : float
42 Initial phase
43 amplitude : float
44 Wave amplitude
45 v_alfven : float
46 Alfvén velocity
47 direction : int
48 +1 for right-propagating, -1 for left-propagating
49
50 Returns
51 -------
52 wave : ndarray
53 Wave packet field
54 """
55 envelope = amplitude * np.exp(-(x - direction * v_alfven * t)**2 / width**2)
56 wave = envelope * np.cos(k0 * (x - direction * v_alfven * t) + phase)
57 return wave
58
59
60def compute_nonlinear_term(vx, vy, Bx, By, dx, v_alfven):
61 """
62 Compute nonlinear term in MHD equations: (v·∇)B - (B·∇)v.
63
64 For simplicity, we use the Elsässer formulation:
65 ∂z±/∂t + v_A ∂z±/∂x = -(z∓·∇)z±
66
67 Parameters
68 ----------
69 vx, vy : ndarray
70 Velocity components
71 Bx, By : ndarray
72 Magnetic field components (in velocity units)
73 dx : float
74 Grid spacing
75 v_alfven : float
76 Alfvén velocity
77
78 Returns
79 -------
80 nl_term : ndarray
81 Nonlinear interaction term
82 """
83 # Elsässer variables
84 zp_x = vx + Bx
85 zp_y = vy + By
86 zm_x = vx - Bx
87 zm_y = vy - By
88
89 # Gradients (simple centered difference)
90 dzm_x_dx = np.gradient(zm_x, dx)
91 dzp_x_dx = np.gradient(zp_x, dx)
92
93 # Nonlinear term: -(z∓·∇)z±
94 nl_zp = -zm_x * dzp_x_dx
95 nl_zm = -zp_x * dzm_x_dx
96
97 return nl_zp, nl_zm
98
99
100def simulate_alfven_collision(N=512, L=100.0, T_max=50.0, dt=0.05):
101 """
102 Simulate collision of two Alfvén wave packets.
103
104 Parameters
105 ----------
106 N : int
107 Number of grid points
108 L : float
109 Domain length
110 T_max : float
111 Maximum simulation time
112 dt : float
113 Time step
114
115 Returns
116 -------
117 x : ndarray
118 Spatial grid
119 t_array : ndarray
120 Time array
121 energy_history : ndarray
122 Total energy vs time
123 spectrum_history : list
124 Energy spectra at different times
125 """
126 # Setup grid
127 x = np.linspace(0, L, N)
128 dx = L / N
129 k = fftfreq(N, d=dx) * 2 * np.pi
130
131 # Parameters
132 v_alfven = 1.0 # Alfvén velocity
133 rho = 1.0 # Density
134
135 # Initial conditions: two counter-propagating wave packets
136 k0 = 2 * np.pi / 10 # Central wavenumber
137 width = 5.0 # Packet width
138 amplitude = 0.3 # Wave amplitude
139
140 # Initialize fields
141 vx = np.zeros(N)
142 vy = (alfven_wave_packet(x, 0, k0, width, 0, amplitude, v_alfven, +1) +
143 alfven_wave_packet(x, 0, k0, width, np.pi/4, amplitude, v_alfven, -1))
144
145 # In Alfvén waves: δv = ±δB/√(ρμ₀), we use normalized units
146 Bx = np.zeros(N)
147 By_right = alfven_wave_packet(x, 0, k0, width, 0, amplitude, v_alfven, +1)
148 By_left = -alfven_wave_packet(x, 0, k0, width, np.pi/4, amplitude, v_alfven, -1)
149 By = By_right + By_left
150
151 # Time integration
152 n_steps = int(T_max / dt)
153 t_array = np.linspace(0, T_max, n_steps)
154 energy_history = np.zeros(n_steps)
155 spectrum_history = []
156 snapshot_times = [0, T_max/4, T_max/2, 3*T_max/4, T_max]
157
158 for step in range(n_steps):
159 t = t_array[step]
160
161 # Compute energy
162 energy = np.sum(vx**2 + vy**2 + Bx**2 + By**2) * dx / 2
163 energy_history[step] = energy
164
165 # Store snapshots for spectrum
166 if any(np.abs(t - st) < dt/2 for st in snapshot_times):
167 vy_k = fft(vy)
168 spectrum = np.abs(vy_k)**2 / N
169 spectrum_history.append((t, k[:N//2], spectrum[:N//2]))
170
171 # Compute nonlinear terms
172 nl_zp, nl_zm = compute_nonlinear_term(vx, vy, Bx, By, dx, v_alfven)
173
174 # Time advancement (simple forward Euler for demonstration)
175 # ∂z+/∂t = -v_A ∂z+/∂x - (z-·∇)z+
176 # ∂z-/∂t = +v_A ∂z-/∂x - (z+·∇)z-
177
178 zp_x = vx + Bx
179 zp_y = vy + By
180 zm_x = vx - Bx
181 zm_y = vy - By
182
183 # Spectral method for linear terms (advection)
184 zp_y_k = fft(zp_y)
185 zm_y_k = fft(zm_y)
186
187 # Advection in Fourier space
188 zp_y_k = zp_y_k * np.exp(-1j * v_alfven * k * dt)
189 zm_y_k = zm_y_k * np.exp(+1j * v_alfven * k * dt)
190
191 # Back to real space
192 zp_y = ifft(zp_y_k).real
193 zm_y = ifft(zm_y_k).real
194
195 # Add nonlinear term
196 zp_y = zp_y + nl_zp * dt
197 zm_y = zm_y + nl_zm * dt
198
199 # Update physical variables
200 vy = 0.5 * (zp_y + zm_y)
201 By = 0.5 * (zp_y - zm_y)
202
203 # Damping at boundaries (prevent reflections)
204 damping = np.exp(-((x - L/2)**2) / (L/3)**2)
205 vy = vy * damping
206 By = By * damping
207
208 return x, t_array, energy_history, spectrum_history
209
210
211def plot_collision_results(x, t_array, energy_history, spectrum_history):
212 """
213 Plot Alfvén wave collision results.
214
215 Parameters
216 ----------
217 x : ndarray
218 Spatial grid
219 t_array : ndarray
220 Time array
221 energy_history : ndarray
222 Energy vs time
223 spectrum_history : list
224 Energy spectra at different times
225 """
226 fig = plt.figure(figsize=(16, 10))
227
228 # Plot 1: Energy conservation
229 ax1 = plt.subplot(2, 3, 1)
230 ax1.plot(t_array, energy_history / energy_history[0], 'b-', linewidth=2)
231 ax1.axhline(1.0, color='k', linestyle='--', linewidth=1, alpha=0.5)
232 ax1.set_xlabel('Time (τ_A)', fontsize=12)
233 ax1.set_ylabel('Normalized Total Energy', fontsize=12)
234 ax1.set_title('Energy Conservation', fontsize=13, fontweight='bold')
235 ax1.grid(True, alpha=0.3)
236
237 # Plot 2-4: Energy spectra at different times
238 colors = ['blue', 'green', 'orange', 'red', 'purple']
239 ax2 = plt.subplot(2, 3, 2)
240
241 for i, (t, k_vals, spectrum) in enumerate(spectrum_history):
242 # Remove zero frequency and plot
243 k_plot = k_vals[1:len(k_vals)//2]
244 spec_plot = spectrum[1:len(spectrum)//2]
245 ax2.loglog(k_plot, spec_plot, color=colors[i % len(colors)],
246 linewidth=2, alpha=0.7, label=f't = {t:.1f}')
247
248 # Add theoretical predictions
249 k_theory = k_vals[2:20]
250 ax2.loglog(k_theory, 100 * k_theory**(-5/3), 'k--',
251 linewidth=1.5, label=r'$k^{-5/3}$', alpha=0.5)
252 ax2.loglog(k_theory, 50 * k_theory**(-3/2), 'k:',
253 linewidth=1.5, label=r'$k^{-3/2}$', alpha=0.5)
254
255 ax2.set_xlabel('Wavenumber k', fontsize=12)
256 ax2.set_ylabel('Energy E(k)', fontsize=12)
257 ax2.set_title('Energy Spectrum Evolution', fontsize=13, fontweight='bold')
258 ax2.legend(loc='upper right', fontsize=9)
259 ax2.grid(True, alpha=0.3, which='both')
260
261 # Plot 3: Critical balance diagram
262 ax3 = plt.subplot(2, 3, 3)
263 k_range = np.logspace(0, 2, 50)
264
265 # Cascade time: τ_nl ~ k_perp / δv_k ~ k_perp / (E(k) k)^(1/2)
266 # For IK: E(k) ~ k^(-3/2), so τ_nl ~ k^(-1/4)
267 tau_nl_ik = 10 * k_range**(-1/4)
268 # Alfvén time: τ_A ~ 1/k_parallel
269 tau_alfven = 10 / k_range
270
271 ax3.loglog(k_range, tau_nl_ik, 'b-', linewidth=2,
272 label=r'$\tau_{nl}$ (cascade time)')
273 ax3.loglog(k_range, tau_alfven, 'r-', linewidth=2,
274 label=r'$\tau_A$ (Alfvén time)')
275 ax3.axhline(10, color='k', linestyle='--', linewidth=1,
276 alpha=0.5, label='Critical balance')
277
278 ax3.set_xlabel('Wavenumber k', fontsize=12)
279 ax3.set_ylabel('Timescale', fontsize=12)
280 ax3.set_title('Critical Balance: τ_nl ~ τ_A', fontsize=13, fontweight='bold')
281 ax3.legend(loc='upper right', fontsize=10)
282 ax3.grid(True, alpha=0.3, which='both')
283
284 # Plot 4: Cascade rate
285 ax4 = plt.subplot(2, 3, 4)
286 # Energy cascade rate: ε ~ E / τ_nl
287 # For critical balance with IK spectrum: ε ~ const
288 k_cascade = k_range[k_range > 1]
289 E_k_kolm = k_cascade**(-5/3)
290 E_k_ik = k_cascade**(-3/2)
291 epsilon_kolm = E_k_kolm / (k_cascade**(-2/3)) # τ_nl ~ k^(-2/3)
292 epsilon_ik = E_k_ik / (k_cascade**(-1/4)) # τ_nl ~ k^(-1/4)
293
294 ax4.semilogx(k_cascade, epsilon_kolm / epsilon_kolm[0], 'b-',
295 linewidth=2, label='Kolmogorov', alpha=0.7)
296 ax4.semilogx(k_cascade, epsilon_ik / epsilon_ik[0], 'r-',
297 linewidth=2, label='IK (critical balance)', alpha=0.7)
298 ax4.axhline(1.0, color='k', linestyle='--', linewidth=1, alpha=0.5)
299
300 ax4.set_xlabel('Wavenumber k', fontsize=12)
301 ax4.set_ylabel('Normalized Cascade Rate ε', fontsize=12)
302 ax4.set_title('Energy Cascade Rate', fontsize=13, fontweight='bold')
303 ax4.legend(loc='best', fontsize=10)
304 ax4.grid(True, alpha=0.3)
305
306 # Plot 5: Anisotropy
307 ax5 = plt.subplot(2, 3, 5)
308 # k_parallel / k_perp ratio for critical balance
309 # For GS95: k_|| ~ k_perp^(2/3)
310 k_perp = np.logspace(0, 2, 50)
311 k_par_gs = k_perp**(2/3) # Goldreich-Sridhar
312 k_par_iso = k_perp # Isotropic
313
314 ax5.loglog(k_perp, k_par_gs, 'b-', linewidth=2,
315 label='GS95: k_|| ~ k_⊥^(2/3)', alpha=0.7)
316 ax5.loglog(k_perp, k_par_iso, 'k--', linewidth=1.5,
317 label='Isotropic: k_|| ~ k_⊥', alpha=0.5)
318
319 ax5.set_xlabel('k_⊥ (perpendicular)', fontsize=12)
320 ax5.set_ylabel('k_|| (parallel)', fontsize=12)
321 ax5.set_title('Spectral Anisotropy', fontsize=13, fontweight='bold')
322 ax5.legend(loc='upper left', fontsize=10)
323 ax5.grid(True, alpha=0.3, which='both')
324
325 # Plot 6: Physics summary
326 ax6 = plt.subplot(2, 3, 6)
327 ax6.axis('off')
328
329 summary_text = """
330 Critical Balance Theory
331 ───────────────────────
332
333 • Cascade mediated by Alfvén wave collisions
334 • Energy transfer occurs when:
335 τ_nl (nonlinear time) ~ τ_A (Alfvén time)
336
337 Spectral Predictions:
338 • Kolmogorov (isotropic): E(k) ~ k^(-5/3)
339 • Iroshnikov-Kraichnan: E(k) ~ k^(-3/2)
340 • Goldreich-Sridhar (anisotropic):
341 - Perpendicular cascade dominates
342 - k_|| ~ k_⊥^(2/3)
343 - E(k_⊥) ~ k_⊥^(-5/3)
344
345 Key Physics:
346 • Counter-propagating Alfvén waves interact
347 • Generate small-scale perturbations
348 • Energy cascades to dissipation scales
349 • Anisotropy due to mean B field
350 """
351
352 ax6.text(0.1, 0.5, summary_text, fontsize=11, family='monospace',
353 verticalalignment='center', transform=ax6.transAxes)
354
355 plt.suptitle('Alfvén Wave Collision and Energy Cascade',
356 fontsize=15, fontweight='bold')
357 plt.tight_layout()
358 plt.savefig('alfven_wave_cascade.png', dpi=150, bbox_inches='tight')
359 plt.show()
360
361
362def main():
363 """Main execution function."""
364 print("="*60)
365 print("Alfvén Wave Collision and Energy Cascade")
366 print("="*60)
367
368 print("\nSimulating counter-propagating Alfvén wave packets...")
369 x, t_array, energy_history, spectrum_history = simulate_alfven_collision(
370 N=512, L=100.0, T_max=40.0, dt=0.05
371 )
372
373 print(f"Simulation complete: {len(t_array)} timesteps")
374 print(f"Energy conservation: {100*(1 - energy_history[-1]/energy_history[0]):.2f}% loss")
375
376 print("\nGenerating plots...")
377 plot_collision_results(x, t_array, energy_history, spectrum_history)
378
379 print("\nPlot saved as 'alfven_wave_cascade.png'")
380
381 print("\n" + "="*60)
382 print("Key Concepts")
383 print("="*60)
384 print("\n1. Critical Balance:")
385 print(" τ_nl ~ τ_A => Energy cascade rate determined by Alfvén waves")
386 print("\n2. Spectral Index:")
387 print(" E(k) ~ k^(-3/2) for strong MHD turbulence (IK theory)")
388 print(" E(k) ~ k^(-5/3) for hydrodynamic turbulence (Kolmogorov)")
389 print("\n3. Anisotropy:")
390 print(" k_|| ~ k_⊥^(2/3) for Goldreich-Sridhar cascade")
391 print("\nReferences:")
392 print(" - Iroshnikov (1964), Kraichnan (1965)")
393 print(" - Goldreich & Sridhar (1995)")
394 print(" - Boldyrev (2006)")
395
396
397if __name__ == '__main__':
398 main()