1#!/usr/bin/env python3
2"""
3MHD Turbulence Energy Spectra
4
5This script computes and visualizes energy spectra in MHD turbulence,
6comparing different theoretical predictions:
7- Kolmogorov k^(-5/3) for hydrodynamic turbulence
8- Iroshnikov-Kraichnan k^(-3/2) for strong MHD turbulence
9- Goldreich-Sridhar for weak turbulence
10
11Key results:
12- Kinetic and magnetic energy spectra show different scalings
13- ElsÀsser variables reveal wave-like nature of MHD turbulence
14- Spectral anisotropy in the presence of a mean magnetic field
15
16Author: Claude
17Date: 2026-02-14
18"""
19
20import numpy as np
21import matplotlib.pyplot as plt
22from scipy.fft import fftn, fftfreq
23from scipy.ndimage import gaussian_filter
24
25
26def generate_turbulent_field(N, alpha, energy_0=1.0, kmin=1, kmax=None):
27 """
28 Generate a turbulent field with power-law spectrum E(k) ~ k^(-alpha).
29
30 Parameters
31 ----------
32 N : int
33 Grid size (cubic grid N^3)
34 alpha : float
35 Spectral index (e.g., 5/3 for Kolmogorov, 3/2 for IK)
36 energy_0 : float
37 Energy normalization
38 kmin : int
39 Minimum wavenumber
40 kmax : int or None
41 Maximum wavenumber (None = N/2)
42
43 Returns
44 -------
45 field : ndarray
46 3D turbulent field
47 """
48 if kmax is None:
49 kmax = N // 2
50
51 # Create 3D field in Fourier space
52 field_k = np.zeros((N, N, N), dtype=complex)
53
54 # Create wavenumber grid
55 kx = fftfreq(N, d=1.0) * N
56 ky = fftfreq(N, d=1.0) * N
57 kz = fftfreq(N, d=1.0) * N
58 KX, KY, KZ = np.meshgrid(kx, ky, kz, indexing='ij')
59 K = np.sqrt(KX**2 + KY**2 + KZ**2)
60
61 # Generate random phases
62 np.random.seed(42)
63 phases = np.random.uniform(0, 2*np.pi, (N, N, N))
64
65 # Set amplitudes according to power law
66 for i in range(N):
67 for j in range(N):
68 for k in range(N):
69 kmag = K[i, j, k]
70 if kmin <= kmag <= kmax and kmag > 0:
71 amplitude = energy_0 * kmag**(-alpha/2)
72 field_k[i, j, k] = amplitude * np.exp(1j * phases[i, j, k])
73
74 # Transform to real space
75 field = np.fft.ifftn(field_k).real
76
77 return field
78
79
80def compute_energy_spectrum_3d(field):
81 """
82 Compute energy spectrum E(k) from a 3D field.
83
84 Parameters
85 ----------
86 field : ndarray
87 3D field (velocity or magnetic)
88
89 Returns
90 -------
91 k_bins : ndarray
92 Wavenumber bins
93 E_k : ndarray
94 Energy spectrum
95 """
96 N = field.shape[0]
97
98 # Fourier transform
99 field_k = fftn(field)
100
101 # Energy in Fourier space: |field_k|^2 / 2
102 energy_k = 0.5 * np.abs(field_k)**2 / N**6
103
104 # Create wavenumber grid
105 kx = fftfreq(N, d=1.0) * N
106 ky = fftfreq(N, d=1.0) * N
107 kz = fftfreq(N, d=1.0) * N
108 KX, KY, KZ = np.meshgrid(kx, ky, kz, indexing='ij')
109 K = np.sqrt(KX**2 + KY**2 + KZ**2)
110
111 # Bin by wavenumber magnitude
112 k_bins = np.arange(1, N//2)
113 E_k = np.zeros(len(k_bins))
114
115 for i, k in enumerate(k_bins):
116 mask = (K >= k - 0.5) & (K < k + 0.5)
117 E_k[i] = np.sum(energy_k[mask])
118
119 return k_bins, E_k
120
121
122def compute_elsasser_variables(vx, vy, vz, Bx, By, Bz, rho=1.0, mu0=1.0):
123 """
124 Compute ElsĂ€sser variables z± = v ± B/â(ÏÎŒâ).
125
126 Parameters
127 ----------
128 vx, vy, vz : ndarray
129 Velocity components
130 Bx, By, Bz : ndarray
131 Magnetic field components
132 rho : float
133 Density (assumed uniform)
134 mu0 : float
135 Magnetic permeability
136
137 Returns
138 -------
139 zp, zm : tuple of ndarray
140 ElsÀsser variables z+ and z-
141 """
142 va_factor = 1.0 / np.sqrt(rho * mu0)
143
144 zp_x = vx + Bx * va_factor
145 zp_y = vy + By * va_factor
146 zp_z = vz + Bz * va_factor
147
148 zm_x = vx - Bx * va_factor
149 zm_y = vy - By * va_factor
150 zm_z = vz - Bz * va_factor
151
152 return (zp_x, zp_y, zp_z), (zm_x, zm_y, zm_z)
153
154
155def plot_energy_spectra():
156 """
157 Plot and compare MHD turbulence energy spectra.
158 """
159 # Grid size
160 N = 64
161
162 # Generate turbulent fields with different spectral indices
163 print("Generating turbulent fields...")
164
165 # Kolmogorov (hydrodynamic)
166 vx_kolm = generate_turbulent_field(N, 5/3, energy_0=1.0)
167 vy_kolm = generate_turbulent_field(N, 5/3, energy_0=1.0)
168 vz_kolm = generate_turbulent_field(N, 5/3, energy_0=1.0)
169
170 # Iroshnikov-Kraichnan (MHD)
171 vx_ik = generate_turbulent_field(N, 3/2, energy_0=1.0)
172 vy_ik = generate_turbulent_field(N, 3/2, energy_0=1.0)
173 vz_ik = generate_turbulent_field(N, 3/2, energy_0=1.0)
174
175 Bx_ik = generate_turbulent_field(N, 3/2, energy_0=0.8)
176 By_ik = generate_turbulent_field(N, 3/2, energy_0=0.8)
177 Bz_ik = generate_turbulent_field(N, 3/2, energy_0=0.8)
178
179 # Compute spectra
180 print("Computing energy spectra...")
181
182 # Kinetic energy spectrum (Kolmogorov)
183 k_bins, E_kin_kolm = compute_energy_spectrum_3d(vx_kolm)
184 _, E_kin_kolm_y = compute_energy_spectrum_3d(vy_kolm)
185 _, E_kin_kolm_z = compute_energy_spectrum_3d(vz_kolm)
186 E_kin_kolm_total = E_kin_kolm + E_kin_kolm_y + E_kin_kolm_z
187
188 # Kinetic and magnetic energy spectra (IK)
189 k_bins, E_kin_ik = compute_energy_spectrum_3d(vx_ik)
190 _, E_kin_ik_y = compute_energy_spectrum_3d(vy_ik)
191 _, E_kin_ik_z = compute_energy_spectrum_3d(vz_ik)
192 E_kin_ik_total = E_kin_ik + E_kin_ik_y + E_kin_ik_z
193
194 _, E_mag_ik = compute_energy_spectrum_3d(Bx_ik)
195 _, E_mag_ik_y = compute_energy_spectrum_3d(By_ik)
196 _, E_mag_ik_z = compute_energy_spectrum_3d(Bz_ik)
197 E_mag_ik_total = E_mag_ik + E_mag_ik_y + E_mag_ik_z
198
199 # ElsÀsser variables
200 zp, zm = compute_elsasser_variables(vx_ik, vy_ik, vz_ik,
201 Bx_ik, By_ik, Bz_ik)
202 _, E_zp = compute_energy_spectrum_3d(zp[0])
203 _, E_zp_y = compute_energy_spectrum_3d(zp[1])
204 _, E_zp_z = compute_energy_spectrum_3d(zp[2])
205 E_zp_total = E_zp + E_zp_y + E_zp_z
206
207 _, E_zm = compute_energy_spectrum_3d(zm[0])
208 _, E_zm_y = compute_energy_spectrum_3d(zm[1])
209 _, E_zm_z = compute_energy_spectrum_3d(zm[2])
210 E_zm_total = E_zm + E_zm_y + E_zm_z
211
212 # Theoretical predictions
213 k_theory = k_bins[k_bins > 2]
214 E_kolm_theory = 100 * k_theory**(-5/3)
215 E_ik_theory = 100 * k_theory**(-3/2)
216
217 # Create plots
218 fig, axes = plt.subplots(2, 2, figsize=(14, 12))
219
220 # Plot 1: Kolmogorov vs Iroshnikov-Kraichnan
221 ax = axes[0, 0]
222 ax.loglog(k_bins, E_kin_kolm_total, 'b-', linewidth=2,
223 label='Kinetic (Kolmogorov)', alpha=0.7)
224 ax.loglog(k_bins, E_kin_ik_total, 'r-', linewidth=2,
225 label='Kinetic (IK)', alpha=0.7)
226 ax.loglog(k_theory, E_kolm_theory, 'b--', linewidth=1.5,
227 label=r'$k^{-5/3}$ (Kolmogorov)')
228 ax.loglog(k_theory, E_ik_theory, 'r--', linewidth=1.5,
229 label=r'$k^{-3/2}$ (IK)')
230 ax.set_xlabel('Wavenumber k', fontsize=12)
231 ax.set_ylabel('Energy E(k)', fontsize=12)
232 ax.set_title('Spectral Index Comparison', fontsize=13, fontweight='bold')
233 ax.legend(loc='upper right', fontsize=10)
234 ax.grid(True, alpha=0.3, which='both')
235
236 # Plot 2: Kinetic vs Magnetic energy
237 ax = axes[0, 1]
238 ax.loglog(k_bins, E_kin_ik_total, 'b-', linewidth=2,
239 label='Kinetic energy', alpha=0.7)
240 ax.loglog(k_bins, E_mag_ik_total, 'r-', linewidth=2,
241 label='Magnetic energy', alpha=0.7)
242 ax.loglog(k_bins, E_kin_ik_total + E_mag_ik_total, 'k-', linewidth=2,
243 label='Total energy', alpha=0.7)
244 ax.loglog(k_theory, E_ik_theory, 'g--', linewidth=1.5,
245 label=r'$k^{-3/2}$')
246 ax.set_xlabel('Wavenumber k', fontsize=12)
247 ax.set_ylabel('Energy E(k)', fontsize=12)
248 ax.set_title('Kinetic vs Magnetic Energy', fontsize=13, fontweight='bold')
249 ax.legend(loc='upper right', fontsize=10)
250 ax.grid(True, alpha=0.3, which='both')
251
252 # Plot 3: ElsÀsser variables
253 ax = axes[1, 0]
254 ax.loglog(k_bins, E_zp_total, 'b-', linewidth=2,
255 label=r'$z^+$ (outgoing)', alpha=0.7)
256 ax.loglog(k_bins, E_zm_total, 'r-', linewidth=2,
257 label=r'$z^-$ (incoming)', alpha=0.7)
258 ax.loglog(k_theory, E_ik_theory, 'k--', linewidth=1.5,
259 label=r'$k^{-3/2}$')
260 ax.set_xlabel('Wavenumber k', fontsize=12)
261 ax.set_ylabel('Energy E(k)', fontsize=12)
262 ax.set_title('ElsĂ€sser Variables z± = v ± B/â(ÏÎŒâ)',
263 fontsize=13, fontweight='bold')
264 ax.legend(loc='upper right', fontsize=10)
265 ax.grid(True, alpha=0.3, which='both')
266
267 # Plot 4: Energy ratio
268 ax = axes[1, 1]
269 ratio = E_mag_ik_total / (E_kin_ik_total + 1e-10)
270 ax.semilogx(k_bins, ratio, 'g-', linewidth=2, alpha=0.7)
271 ax.axhline(1.0, color='k', linestyle='--', linewidth=1.5,
272 label='Equipartition')
273 ax.set_xlabel('Wavenumber k', fontsize=12)
274 ax.set_ylabel('E_mag / E_kin', fontsize=12)
275 ax.set_title('Energy Equipartition', fontsize=13, fontweight='bold')
276 ax.set_ylim([0, 3])
277 ax.legend(loc='upper right', fontsize=10)
278 ax.grid(True, alpha=0.3)
279
280 plt.suptitle('MHD Turbulence Energy Spectra',
281 fontsize=15, fontweight='bold')
282 plt.tight_layout()
283 plt.savefig('mhd_energy_spectrum.png', dpi=150, bbox_inches='tight')
284 plt.show()
285
286 # Print summary
287 print("\n" + "="*60)
288 print("MHD Turbulence Energy Spectra")
289 print("="*60)
290 print("\nSpectral Indices:")
291 print(" Kolmogorov (hydrodynamic): E(k) ~ k^(-5/3) â k^(-1.67)")
292 print(" Iroshnikov-Kraichnan (MHD): E(k) ~ k^(-3/2) = k^(-1.50)")
293 print("\nElsÀsser Variables:")
294 print(" z+ = v + B/â(ÏÎŒâ) (outgoing AlfvĂ©n waves)")
295 print(" z- = v - B/â(ÏÎŒâ) (incoming AlfvĂ©n waves)")
296 print("\nKey Physics:")
297 print(" - MHD turbulence is mediated by Alfvén wave collisions")
298 print(" - IK theory predicts shallower spectrum than Kolmogorov")
299 print(" - Energy equipartition: E_kin â E_mag in inertial range")
300 print(" - Critical balance: energy cascade rate ~ Alfvén wave period")
301
302
303def main():
304 """Main execution function."""
305 print("MHD Turbulence Energy Spectrum Analysis")
306 print("=" * 60)
307
308 plot_energy_spectra()
309
310 print("\nPlot saved as 'mhd_energy_spectrum.png'")
311 print("\nReferences:")
312 print(" - Kolmogorov (1941): k^(-5/3) for hydrodynamic turbulence")
313 print(" - Iroshnikov (1964), Kraichnan (1965): k^(-3/2) for MHD")
314 print(" - Goldreich & Sridhar (1995): Anisotropic cascade theory")
315
316
317if __name__ == '__main__':
318 main()