1#!/usr/bin/env python3
2"""
3Structure Functions and Intermittency in MHD Turbulence
4
5This script computes structure functions of velocity and magnetic field
6to characterize the scaling properties and intermittency of MHD turbulence.
7
8Structure functions:
9 S_p(ℓ) = <|δv(ℓ)|^p>
10 B_p(ℓ) = <|δB(ℓ)|^p>
11
12where δv(ℓ) = v(x + ℓ) - v(x) is the velocity increment.
13
14Key results:
15- Scaling exponents ζ(p) deviate from Kolmogorov prediction p/3
16- Intermittency: rare, intense fluctuations at small scales
17- Probability density functions (PDFs) show non-Gaussian tails
18- Magnetic field shows similar intermittency to velocity
19
20Author: Claude
21Date: 2026-02-14
22"""
23
24import numpy as np
25import matplotlib.pyplot as plt
26from scipy.special import gamma
27from scipy.optimize import curve_fit
28from scipy.stats import kurtosis, skew
29
30
31def generate_multifractal_field(N, alpha, intermittency=0.3):
32 """
33 Generate a multifractal field with intermittency.
34
35 Uses a multiplicative cascade model to create fields with
36 non-Gaussian statistics and scale-dependent structure.
37
38 Parameters
39 ----------
40 N : int
41 Number of grid points
42 alpha : float
43 Spectral index
44 intermittency : float
45 Intermittency parameter (0 = no intermittency)
46
47 Returns
48 -------
49 field : ndarray
50 1D field with intermittent structures
51 """
52 # Generate Gaussian random field with power-law spectrum
53 np.random.seed(42)
54 k = np.fft.fftfreq(N, d=1.0) * N
55 k[0] = 1e-10 # Avoid division by zero
56
57 # Power spectrum
58 P_k = np.abs(k)**(-alpha)
59 P_k[0] = 0
60
61 # Random phases
62 phases = np.random.uniform(0, 2*np.pi, N)
63 field_k = np.sqrt(P_k) * np.exp(1j * phases)
64
65 # Transform to real space
66 field = np.fft.ifft(field_k).real
67
68 # Add intermittency via multiplicative cascade
69 if intermittency > 0:
70 # Generate random multipliers
71 n_levels = int(np.log2(N))
72 multipliers = np.ones(N)
73
74 for level in range(n_levels):
75 block_size = 2**(n_levels - level)
76 n_blocks = N // block_size
77
78 for i in range(n_blocks):
79 # Log-normal multiplier
80 sigma = intermittency / np.sqrt(level + 1)
81 mult = np.exp(np.random.normal(0, sigma) - sigma**2/2)
82 multipliers[i*block_size:(i+1)*block_size] *= mult
83
84 field = field * multipliers
85
86 # Normalize
87 field = field / np.std(field)
88
89 return field
90
91
92def compute_structure_function(field, order, max_lag=None):
93 """
94 Compute structure function S_p(ℓ) = <|δfield(ℓ)|^p>.
95
96 Parameters
97 ----------
98 field : ndarray
99 1D field (velocity or magnetic)
100 order : float
101 Order p of structure function
102 max_lag : int or None
103 Maximum lag (None = N/4)
104
105 Returns
106 -------
107 lags : ndarray
108 Lag distances ℓ
109 S_p : ndarray
110 Structure function values
111 """
112 N = len(field)
113 if max_lag is None:
114 max_lag = N // 4
115
116 lags = np.arange(1, max_lag)
117 S_p = np.zeros(len(lags))
118
119 for i, lag in enumerate(lags):
120 # Compute increments
121 delta_field = field[lag:] - field[:-lag]
122 # Structure function: average of |δfield|^p
123 S_p[i] = np.mean(np.abs(delta_field)**order)
124
125 return lags, S_p
126
127
128def compute_scaling_exponent(lags, S_p, fit_range=(0.1, 0.5)):
129 """
130 Compute scaling exponent ζ(p) from S_p(ℓ) ~ ℓ^ζ(p).
131
132 Parameters
133 ----------
134 lags : ndarray
135 Lag distances
136 S_p : ndarray
137 Structure function
138 fit_range : tuple
139 (min_fraction, max_fraction) of lags to use for fitting
140
141 Returns
142 -------
143 zeta : float
144 Scaling exponent
145 """
146 # Select fitting range
147 N = len(lags)
148 i_min = int(fit_range[0] * N)
149 i_max = int(fit_range[1] * N)
150
151 # Log-log fit
152 log_lags = np.log(lags[i_min:i_max])
153 log_S_p = np.log(S_p[i_min:i_max])
154
155 # Linear fit
156 coeffs = np.polyfit(log_lags, log_S_p, 1)
157 zeta = coeffs[0]
158
159 return zeta
160
161
162def compute_pdf(field, n_bins=50):
163 """
164 Compute probability density function of field increments.
165
166 Parameters
167 ----------
168 field : ndarray
169 1D field
170 n_bins : int
171 Number of bins for histogram
172
173 Returns
174 -------
175 bins : ndarray
176 Bin centers
177 pdf : ndarray
178 Probability density
179 """
180 # Compute increments at small scale
181 lag = len(field) // 20
182 delta_field = field[lag:] - field[:-lag]
183
184 # Normalize by standard deviation
185 delta_field = delta_field / np.std(delta_field)
186
187 # Compute histogram
188 counts, bin_edges = np.histogram(delta_field, bins=n_bins, density=True)
189 bins = 0.5 * (bin_edges[1:] + bin_edges[:-1])
190
191 return bins, counts
192
193
194def gaussian_pdf(x):
195 """Standard Gaussian PDF."""
196 return np.exp(-x**2 / 2) / np.sqrt(2 * np.pi)
197
198
199def plot_structure_functions():
200 """
201 Compute and plot structure functions for MHD turbulence.
202 """
203 # Generate turbulent fields
204 N = 2048
205 print("Generating turbulent fields...")
206
207 # Velocity field (with intermittency)
208 v_field = generate_multifractal_field(N, alpha=5/3, intermittency=0.4)
209
210 # Magnetic field (with intermittency)
211 B_field = generate_multifractal_field(N, alpha=5/3, intermittency=0.4)
212
213 # Gaussian field (no intermittency, for comparison)
214 v_gaussian = generate_multifractal_field(N, alpha=5/3, intermittency=0.0)
215
216 # Compute structure functions for different orders
217 orders = [1, 2, 3, 4, 5, 6]
218 print("Computing structure functions...")
219
220 fig = plt.figure(figsize=(16, 12))
221
222 # Plot 1: Structure functions for velocity
223 ax1 = plt.subplot(3, 3, 1)
224 colors = plt.cm.viridis(np.linspace(0, 1, len(orders)))
225
226 for i, p in enumerate(orders):
227 lags, S_p = compute_structure_function(v_field, p)
228 ax1.loglog(lags, S_p, 'o-', color=colors[i], markersize=3,
229 linewidth=1.5, alpha=0.7, label=f'p={p}')
230
231 # Kolmogorov prediction: S_p ~ ℓ^(p/3)
232 if p <= 3:
233 lags_theory = lags[::10]
234 S_theory = lags_theory**(p/3) * S_p[5] / lags[5]**(p/3)
235 ax1.loglog(lags_theory, S_theory, '--', color=colors[i],
236 linewidth=1, alpha=0.5)
237
238 ax1.set_xlabel('Lag ℓ', fontsize=12)
239 ax1.set_ylabel('S_p(ℓ)', fontsize=12)
240 ax1.set_title('Velocity Structure Functions', fontsize=13, fontweight='bold')
241 ax1.legend(loc='lower right', fontsize=9)
242 ax1.grid(True, alpha=0.3, which='both')
243
244 # Plot 2: Structure functions for magnetic field
245 ax2 = plt.subplot(3, 3, 2)
246
247 for i, p in enumerate(orders):
248 lags, B_p = compute_structure_function(B_field, p)
249 ax2.loglog(lags, B_p, 'o-', color=colors[i], markersize=3,
250 linewidth=1.5, alpha=0.7, label=f'p={p}')
251
252 ax2.set_xlabel('Lag ℓ', fontsize=12)
253 ax2.set_ylabel('B_p(ℓ)', fontsize=12)
254 ax2.set_title('Magnetic Structure Functions', fontsize=13, fontweight='bold')
255 ax2.legend(loc='lower right', fontsize=9)
256 ax2.grid(True, alpha=0.3, which='both')
257
258 # Plot 3: Scaling exponents ζ(p)
259 ax3 = plt.subplot(3, 3, 3)
260 orders_full = np.arange(1, 9)
261 zeta_v = np.zeros(len(orders_full))
262 zeta_B = np.zeros(len(orders_full))
263 zeta_gaussian = np.zeros(len(orders_full))
264
265 for i, p in enumerate(orders_full):
266 lags, S_p_v = compute_structure_function(v_field, p)
267 lags, S_p_B = compute_structure_function(B_field, p)
268 lags, S_p_g = compute_structure_function(v_gaussian, p)
269
270 zeta_v[i] = compute_scaling_exponent(lags, S_p_v)
271 zeta_B[i] = compute_scaling_exponent(lags, S_p_B)
272 zeta_gaussian[i] = compute_scaling_exponent(lags, S_p_g)
273
274 # Kolmogorov prediction: ζ(p) = p/3
275 zeta_kolm = orders_full / 3
276
277 ax3.plot(orders_full, zeta_v, 'bo-', linewidth=2, markersize=6,
278 label='Velocity (intermittent)', alpha=0.7)
279 ax3.plot(orders_full, zeta_B, 'rs-', linewidth=2, markersize=6,
280 label='Magnetic (intermittent)', alpha=0.7)
281 ax3.plot(orders_full, zeta_gaussian, 'g^-', linewidth=2, markersize=6,
282 label='Gaussian (no intermittency)', alpha=0.7)
283 ax3.plot(orders_full, zeta_kolm, 'k--', linewidth=2,
284 label='Kolmogorov: ζ = p/3', alpha=0.5)
285
286 ax3.set_xlabel('Order p', fontsize=12)
287 ax3.set_ylabel('Scaling Exponent ζ(p)', fontsize=12)
288 ax3.set_title('Scaling Exponents vs Order', fontsize=13, fontweight='bold')
289 ax3.legend(loc='upper left', fontsize=10)
290 ax3.grid(True, alpha=0.3)
291
292 # Plot 4: Intermittency measure
293 ax4 = plt.subplot(3, 3, 4)
294 # Intermittency: deviation from Kolmogorov
295 Delta_zeta_v = zeta_kolm - zeta_v
296 Delta_zeta_B = zeta_kolm - zeta_B
297
298 ax4.plot(orders_full, Delta_zeta_v, 'bo-', linewidth=2, markersize=6,
299 label='Velocity', alpha=0.7)
300 ax4.plot(orders_full, Delta_zeta_B, 'rs-', linewidth=2, markersize=6,
301 label='Magnetic', alpha=0.7)
302 ax4.axhline(0, color='k', linestyle='--', linewidth=1, alpha=0.5)
303
304 ax4.set_xlabel('Order p', fontsize=12)
305 ax4.set_ylabel('Δζ(p) = p/3 - ζ(p)', fontsize=12)
306 ax4.set_title('Intermittency Correction', fontsize=13, fontweight='bold')
307 ax4.legend(loc='upper left', fontsize=10)
308 ax4.grid(True, alpha=0.3)
309
310 # Plot 5: PDF of velocity increments
311 ax5 = plt.subplot(3, 3, 5)
312
313 x_bins_v, pdf_v = compute_pdf(v_field)
314 x_bins_B, pdf_B = compute_pdf(B_field)
315 x_bins_g, pdf_g = compute_pdf(v_gaussian)
316
317 x_theory = np.linspace(-5, 5, 100)
318 pdf_theory = gaussian_pdf(x_theory)
319
320 ax5.semilogy(x_bins_v, pdf_v, 'b-', linewidth=2,
321 label='Velocity (intermittent)', alpha=0.7)
322 ax5.semilogy(x_bins_g, pdf_g, 'g-', linewidth=2,
323 label='Gaussian field', alpha=0.7)
324 ax5.semilogy(x_theory, pdf_theory, 'k--', linewidth=2,
325 label='Gaussian', alpha=0.5)
326
327 ax5.set_xlabel('δv / σ', fontsize=12)
328 ax5.set_ylabel('PDF', fontsize=12)
329 ax5.set_title('Velocity Increment PDF', fontsize=13, fontweight='bold')
330 ax5.set_ylim([1e-4, 10])
331 ax5.legend(loc='upper right', fontsize=10)
332 ax5.grid(True, alpha=0.3)
333
334 # Plot 6: PDF of magnetic increments
335 ax6 = plt.subplot(3, 3, 6)
336
337 ax6.semilogy(x_bins_B, pdf_B, 'r-', linewidth=2,
338 label='Magnetic (intermittent)', alpha=0.7)
339 ax6.semilogy(x_theory, pdf_theory, 'k--', linewidth=2,
340 label='Gaussian', alpha=0.5)
341
342 ax6.set_xlabel('δB / σ', fontsize=12)
343 ax6.set_ylabel('PDF', fontsize=12)
344 ax6.set_title('Magnetic Increment PDF', fontsize=13, fontweight='bold')
345 ax6.set_ylim([1e-4, 10])
346 ax6.legend(loc='upper right', fontsize=10)
347 ax6.grid(True, alpha=0.3)
348
349 # Plot 7: Higher-order statistics
350 ax7 = plt.subplot(3, 3, 7)
351
352 # Compute statistics at different scales
353 n_scales = 20
354 scales = np.logspace(0, np.log10(N//4), n_scales).astype(int)
355 kurt_v = np.zeros(n_scales)
356 kurt_B = np.zeros(n_scales)
357 skew_v = np.zeros(n_scales)
358
359 for i, scale in enumerate(scales):
360 delta_v = v_field[scale:] - v_field[:-scale]
361 delta_B = B_field[scale:] - B_field[:-scale]
362 kurt_v[i] = kurtosis(delta_v)
363 kurt_B[i] = kurtosis(delta_B)
364 skew_v[i] = skew(delta_v)
365
366 ax7.semilogx(scales, kurt_v, 'b-', linewidth=2, markersize=5,
367 label='Kurtosis (velocity)', alpha=0.7)
368 ax7.semilogx(scales, kurt_B, 'r-', linewidth=2, markersize=5,
369 label='Kurtosis (magnetic)', alpha=0.7)
370 ax7.axhline(0, color='k', linestyle='--', linewidth=1, alpha=0.5,
371 label='Gaussian value')
372
373 ax7.set_xlabel('Scale ℓ', fontsize=12)
374 ax7.set_ylabel('Kurtosis', fontsize=12)
375 ax7.set_title('Scale-Dependent Kurtosis', fontsize=13, fontweight='bold')
376 ax7.legend(loc='upper right', fontsize=10)
377 ax7.grid(True, alpha=0.3)
378
379 # Plot 8: Skewness
380 ax8 = plt.subplot(3, 3, 8)
381
382 ax8.semilogx(scales, skew_v, 'b-', linewidth=2, alpha=0.7)
383 ax8.axhline(0, color='k', linestyle='--', linewidth=1, alpha=0.5)
384
385 ax8.set_xlabel('Scale ℓ', fontsize=12)
386 ax8.set_ylabel('Skewness', fontsize=12)
387 ax8.set_title('Scale-Dependent Skewness', fontsize=13, fontweight='bold')
388 ax8.grid(True, alpha=0.3)
389
390 # Plot 9: Summary text
391 ax9 = plt.subplot(3, 3, 9)
392 ax9.axis('off')
393
394 summary_text = f"""
395 Intermittency Statistics
396 ────────────────────────
397
398 Velocity Field:
399 • Kurtosis (small scale): {kurt_v[0]:.2f}
400 • Skewness: {skew_v[0]:.2f}
401 • ζ(3) deviation: {Delta_zeta_v[2]:.3f}
402
403 Magnetic Field:
404 • Kurtosis (small scale): {kurt_B[0]:.2f}
405 • ζ(3) deviation: {Delta_zeta_B[2]:.3f}
406
407 Key Observations:
408 • Non-Gaussian PDFs
409 • Scale-dependent statistics
410 • ζ(p) < p/3 for p > 3
411 • Enhanced small-scale fluctuations
412
413 Models:
414 • She-Leveque (1994)
415 • Log-normal cascades
416 • Multifractal formalism
417 """
418
419 ax9.text(0.1, 0.5, summary_text, fontsize=10, family='monospace',
420 verticalalignment='center', transform=ax9.transAxes)
421
422 plt.suptitle('Structure Functions and Intermittency in MHD Turbulence',
423 fontsize=15, fontweight='bold')
424 plt.tight_layout()
425 plt.savefig('structure_functions.png', dpi=150, bbox_inches='tight')
426 plt.show()
427
428 # Print numerical results
429 print("\n" + "="*60)
430 print("Scaling Exponents ζ(p)")
431 print("="*60)
432 print(f"{'Order p':<10} {'Velocity':<12} {'Magnetic':<12} {'Kolmogorov':<12}")
433 print("-"*60)
434 for i, p in enumerate(orders_full):
435 print(f"{p:<10} {zeta_v[i]:<12.4f} {zeta_B[i]:<12.4f} {zeta_kolm[i]:<12.4f}")
436
437 print("\n" + "="*60)
438 print("Higher-Order Statistics (small scales)")
439 print("="*60)
440 print(f"Velocity kurtosis: {kurt_v[0]:.3f} (Gaussian = 0)")
441 print(f"Velocity skewness: {skew_v[0]:.3f} (Gaussian = 0)")
442 print(f"Magnetic kurtosis: {kurt_B[0]:.3f} (Gaussian = 0)")
443 print("\nInterpretation:")
444 print(" Positive kurtosis => Heavy tails (intermittency)")
445 print(" Negative skewness => Asymmetry (downward spikes)")
446
447
448def main():
449 """Main execution function."""
450 print("="*60)
451 print("Structure Functions and Intermittency Analysis")
452 print("="*60)
453
454 plot_structure_functions()
455
456 print("\nPlot saved as 'structure_functions.png'")
457
458 print("\n" + "="*60)
459 print("Theoretical Background")
460 print("="*60)
461 print("\nKolmogorov 1941 (K41):")
462 print(" S_p(ℓ) ~ ℓ^(p/3) (self-similar scaling)")
463 print("\nIntermittency Corrections:")
464 print(" S_p(ℓ) ~ ℓ^ζ(p) with ζ(p) < p/3 for p > 3")
465 print(" Due to spatially concentrated dissipation")
466 print("\nMHD Turbulence:")
467 print(" Similar intermittency in velocity and magnetic fields")
468 print(" Current sheets, magnetic islands enhance small-scale activity")
469 print("\nReferences:")
470 print(" - Kolmogorov (1941, 1962)")
471 print(" - She & Leveque (1994)")
472 print(" - Frisch (1995) - Turbulence textbook")
473
474
475if __name__ == '__main__':
476 main()