1#!/usr/bin/env python3
2"""
3Alpha-Omega Dynamo in Spherical Shell
4
5Mean-field dynamo model for magnetic field generation in rotating spheres
6(stars, planets). Demonstrates how differential rotation (Ω-effect) and
7helical turbulence (α-effect) sustain magnetic fields against Ohmic decay.
8
9Key results:
10- Self-sustained magnetic field oscillations
11- Butterfly diagram showing equatorward migration of magnetic activity
12- Periodic field reversals (like Earth's magnetic field)
13- Critical dynamo number for field generation
14
15Physics:
16- Mean-field induction equation with α and Ω effects
17- ∂B/∂t = ∇×(α B + Ω×r × B) + η ∇²B
18- Spherical shell geometry (inner core + convecting outer layer)
19
20Author: Claude
21Date: 2026-02-14
22"""
23
24import numpy as np
25import matplotlib.pyplot as plt
26from scipy.special import sph_harm
27from scipy.integrate import odeint
28from typing import Tuple, List
29
30
31class AlphaOmegaDynamo:
32 """
33 Alpha-Omega dynamo in a spherical shell.
34
35 Simplified mean-field model:
36 - α-effect: Twisting of toroidal field to poloidal field by helical convection
37 - Ω-effect: Shearing of poloidal field to toroidal field by differential rotation
38
39 We use a spectral method with low-order spherical harmonics.
40 """
41
42 def __init__(self, r_inner: float = 0.3, r_outer: float = 1.0,
43 nr: int = 50, ntheta: int = 40):
44 """
45 Initialize dynamo model.
46
47 Args:
48 r_inner: Inner boundary (solid inner core)
49 r_outer: Outer boundary (surface)
50 nr: Radial grid points
51 ntheta: Latitudinal grid points
52 """
53 self.r_inner = r_inner
54 self.r_outer = r_outer
55 self.nr = nr
56 self.ntheta = ntheta
57
58 # Grid
59 self.r = np.linspace(r_inner, r_outer, nr)
60 self.theta = np.linspace(0, np.pi, ntheta) # Colatitude
61 self.R, self.Theta = np.meshgrid(self.r, self.theta, indexing='ij')
62
63 # Latitude (for plotting)
64 self.lat = 90 - np.degrees(self.theta)
65
66 # Time history
67 self.time_history = []
68 self.Bp_history = [] # Poloidal field
69 self.Bt_history = [] # Toroidal field
70
71 def alpha_profile(self, r: np.ndarray, theta: np.ndarray) -> np.ndarray:
72 """
73 Alpha effect profile: α(r, θ)
74
75 Models helical turbulence from convection.
76 Typically: α ∝ cos(θ) (antisymmetric about equator)
77
78 Args:
79 r: Radial coordinate (normalized)
80 theta: Colatitude
81
82 Returns:
83 α coefficient
84 """
85 # Radial profile: peaked in convection zone
86 r_norm = (r - self.r_inner) / (self.r_outer - self.r_inner)
87 radial = np.sin(np.pi * r_norm) # Peaked in middle of shell
88
89 # Latitudinal profile: cos(θ) (north-south antisymmetry)
90 latitudinal = np.cos(theta)
91
92 return radial * latitudinal
93
94 def omega_profile(self, r: np.ndarray, theta: np.ndarray) -> np.ndarray:
95 """
96 Differential rotation profile: Ω(r, θ)
97
98 Solar-like differential rotation:
99 - Faster rotation at equator than poles
100 - Radial shear
101
102 Args:
103 r: Radial coordinate
104 theta: Colatitude
105
106 Returns:
107 Angular velocity Ω
108 """
109 # Latitudinal differential rotation: Ω = Ω₀(1 - β sin²θ)
110 # where β ~ 0.2 for the Sun
111 beta = 0.2
112 lat_term = 1 - beta * np.sin(theta)**2
113
114 # Radial shear: Ω increases outward in convection zone
115 r_norm = (r - self.r_inner) / (self.r_outer - self.r_inner)
116 rad_term = 1 + 0.3 * r_norm
117
118 return lat_term * rad_term
119
120 def initial_field_perturbation(self, amplitude: float = 0.1) -> Tuple[np.ndarray, np.ndarray]:
121 """
122 Create small seed magnetic field to start the dynamo.
123
124 Returns:
125 (B_poloidal, B_toroidal) on the grid
126 """
127 # Small random dipole-like poloidal field
128 # Using P₁¹(cos θ) = sin θ spherical harmonic
129 Bp = amplitude * self.R * np.sin(self.Theta) * np.random.randn()
130
131 # No initial toroidal field
132 Bt = np.zeros_like(Bp)
133
134 return Bp, Bt
135
136 def rhs_dynamo(self, t: float, state: np.ndarray,
137 C_alpha: float, C_omega: float, eta: float) -> np.ndarray:
138 """
139 Right-hand side of the mean-field dynamo equations.
140
141 Equations (simplified, 1D radial + latitudinal):
142 ∂B_p/∂t = C_α B_t + η ∇²B_p
143 ∂B_t/∂t = C_Ω (∇B_p)·(∇Ω) + η ∇²B_t
144
145 Where:
146 - B_p: Poloidal field component
147 - B_t: Toroidal field component
148 - C_α, C_Ω: Dynamo numbers (control strength of α and Ω effects)
149 - η: Magnetic diffusivity
150
151 Args:
152 t: Time
153 state: Flattened [B_p, B_t]
154 C_alpha: Alpha effect strength
155 C_omega: Omega effect strength
156 eta: Magnetic diffusivity
157
158 Returns:
159 d(state)/dt
160 """
161 n = len(state) // 2
162 Bp = state[:n].reshape(self.nr, self.ntheta)
163 Bt = state[n:].reshape(self.nr, self.ntheta)
164
165 # Alpha and Omega profiles
166 alpha = C_alpha * self.alpha_profile(self.R, self.Theta)
167 omega = C_omega * self.omega_profile(self.R, self.Theta)
168
169 # 1. Poloidal field evolution: ∂B_p/∂t = α B_t + η ∇²B_p
170 # Alpha effect: toroidal → poloidal
171 source_Bp = alpha * Bt
172
173 # Diffusion (simplified radial diffusion only)
174 dr = self.r[1] - self.r[0]
175 diff_Bp = np.zeros_like(Bp)
176 diff_Bp[1:-1, :] = eta * (Bp[2:, :] - 2*Bp[1:-1, :] + Bp[:-2, :]) / dr**2
177
178 dBp_dt = source_Bp + diff_Bp
179
180 # Boundary conditions: B_p = 0 at inner and outer boundaries
181 dBp_dt[0, :] = -Bp[0, :] / 0.01 # Relax to zero
182 dBp_dt[-1, :] = -Bp[-1, :] / 0.01
183
184 # 2. Toroidal field evolution: ∂B_t/∂t = Ω-shear * B_p + η ∇²B_t
185 # Omega effect: differential rotation shears poloidal → toroidal
186 # Simplified: dB_t/dt ∝ Ω * ∂B_p/∂r
187 dBp_dr = np.gradient(Bp, dr, axis=0)
188 source_Bt = omega * dBp_dr
189
190 # Diffusion
191 diff_Bt = np.zeros_like(Bt)
192 diff_Bt[1:-1, :] = eta * (Bt[2:, :] - 2*Bt[1:-1, :] + Bt[:-2, :]) / dr**2
193
194 dBt_dt = source_Bt + diff_Bt
195
196 # Boundary conditions
197 dBt_dt[0, :] = -Bt[0, :] / 0.01
198 dBt_dt[-1, :] = -Bt[-1, :] / 0.01
199
200 # Flatten and return
201 return np.concatenate([dBp_dt.flatten(), dBt_dt.flatten()])
202
203 def run_dynamo(self, t_final: float = 100.0, dt: float = 0.1,
204 C_alpha: float = 1.0, C_omega: float = 5.0,
205 eta: float = 0.1) -> None:
206 """
207 Run the dynamo simulation.
208
209 Args:
210 t_final: Final time
211 dt: Time step
212 C_alpha: Alpha effect strength (dynamo number)
213 C_omega: Omega effect strength
214 eta: Magnetic diffusivity
215 """
216 print(f"Running alpha-omega dynamo...")
217 print(f" Parameters: C_α={C_alpha}, C_Ω={C_omega}, η={eta}")
218 print(f" Dynamo number D = C_α * C_Ω = {C_alpha * C_omega:.1f}")
219
220 # Initial conditions
221 Bp0, Bt0 = self.initial_field_perturbation(amplitude=0.01)
222 state0 = np.concatenate([Bp0.flatten(), Bt0.flatten()])
223
224 # Time points
225 t_eval = np.arange(0, t_final + dt, dt)
226
227 # Integrate using explicit Euler (simple for demonstration)
228 self.time_history = [0]
229 self.Bp_history = [Bp0]
230 self.Bt_history = [Bt0]
231
232 state = state0.copy()
233 for i, t in enumerate(t_eval[1:]):
234 # Euler step
235 dstate = self.rhs_dynamo(t, state, C_alpha, C_omega, eta)
236 state = state + dt * dstate
237
238 # Store
239 if (i + 1) % 10 == 0: # Store every 10 steps
240 n = len(state) // 2
241 Bp = state[:n].reshape(self.nr, self.ntheta)
242 Bt = state[n:].reshape(self.nr, self.ntheta)
243
244 self.time_history.append(t)
245 self.Bp_history.append(Bp.copy())
246 self.Bt_history.append(Bt.copy())
247
248 if (i + 1) % 100 == 0:
249 Bp_max = np.max(np.abs(Bp))
250 Bt_max = np.max(np.abs(Bt))
251 print(f" t={t:.1f}, |Bp|_max={Bp_max:.3e}, |Bt|_max={Bt_max:.3e}")
252
253 print(f"Dynamo simulation complete!")
254
255
256def visualize_dynamo(dynamo: AlphaOmegaDynamo, save_prefix: str = "dynamo"):
257 """
258 Visualize dynamo results: butterfly diagram and field snapshots.
259 """
260 times = np.array(dynamo.time_history)
261 n_times = len(times)
262
263 # =========================================================================
264 # 1. Butterfly Diagram (Hovmöller plot)
265 # =========================================================================
266 fig, axes = plt.subplots(3, 2, figsize=(14, 12))
267
268 # Extract toroidal field at mid-radius as function of time and latitude
269 r_mid_idx = dynamo.nr // 2
270 butterfly_data = np.zeros((n_times, dynamo.ntheta))
271
272 for i, Bt in enumerate(dynamo.Bt_history):
273 butterfly_data[i, :] = Bt[r_mid_idx, :]
274
275 # Butterfly diagram
276 ax = axes[0, 0]
277 lat = 90 - np.degrees(dynamo.theta)
278 extent = [lat.min(), lat.max(), times.min(), times.max()]
279 im = ax.imshow(butterfly_data, aspect='auto', cmap='RdBu_r',
280 extent=extent, origin='lower', interpolation='bilinear')
281 ax.set_xlabel('Latitude (degrees)')
282 ax.set_ylabel('Time')
283 ax.set_title('Butterfly Diagram (Toroidal Field B_t)')
284 ax.axvline(x=0, color='k', linestyle='--', alpha=0.5)
285 plt.colorbar(im, ax=ax, label='B_t')
286
287 # =========================================================================
288 # 2. Time series of field strength
289 # =========================================================================
290 ax = axes[0, 1]
291
292 Bp_max_time = [np.max(np.abs(Bp)) for Bp in dynamo.Bp_history]
293 Bt_max_time = [np.max(np.abs(Bt)) for Bt in dynamo.Bt_history]
294
295 ax.plot(times, Bp_max_time, 'b-', label='|Bp|_max', linewidth=2)
296 ax.plot(times, Bt_max_time, 'r-', label='|Bt|_max', linewidth=2)
297 ax.set_xlabel('Time')
298 ax.set_ylabel('Field strength')
299 ax.set_title('Magnetic Field Evolution')
300 ax.legend()
301 ax.grid(True, alpha=0.3)
302 ax.set_yscale('log')
303
304 # =========================================================================
305 # 3. Snapshots at different times
306 # =========================================================================
307 snapshot_indices = [len(times)//4, len(times)//2, 3*len(times)//4, -1]
308
309 for idx, snap_idx in enumerate(snapshot_indices):
310 if idx >= 4:
311 break
312
313 row = 1 + idx // 2
314 col = idx % 2
315 ax = axes[row, col]
316
317 Bt = dynamo.Bt_history[snap_idx]
318 time_snap = times[snap_idx]
319
320 # Plot toroidal field in meridional plane
321 lat_2d = 90 - np.degrees(dynamo.Theta)
322 levels = np.linspace(-np.abs(Bt).max(), np.abs(Bt).max(), 21)
323 cf = ax.contourf(lat_2d, dynamo.R, Bt, levels=levels, cmap='RdBu_r', extend='both')
324 ax.contour(lat_2d, dynamo.R, Bt, levels=levels, colors='k',
325 linewidths=0.5, alpha=0.3)
326
327 ax.set_xlabel('Latitude (degrees)')
328 ax.set_ylabel('Radius')
329 ax.set_title(f'Toroidal Field B_t at t={time_snap:.1f}')
330 ax.axvline(x=0, color='k', linestyle='--', alpha=0.5)
331 plt.colorbar(cf, ax=ax, label='B_t')
332
333 plt.tight_layout()
334 plt.savefig(f'/opt/projects/01_Personal/03_Study/examples/Numerical_Simulation/10_projects/{save_prefix}_butterfly.png',
335 dpi=150, bbox_inches='tight')
336 plt.close()
337 print(f"Figure saved: {save_prefix}_butterfly.png")
338
339 # =========================================================================
340 # 4. Field Structure at Final Time
341 # =========================================================================
342 fig, axes = plt.subplots(1, 2, figsize=(14, 6))
343
344 Bp_final = dynamo.Bp_history[-1]
345 Bt_final = dynamo.Bt_history[-1]
346 lat_2d = 90 - np.degrees(dynamo.Theta)
347
348 # Poloidal field
349 ax = axes[0]
350 levels_p = np.linspace(-np.abs(Bp_final).max(), np.abs(Bp_final).max(), 21)
351 cf = ax.contourf(lat_2d, dynamo.R, Bp_final, levels=levels_p, cmap='PRGn', extend='both')
352 ax.contour(lat_2d, dynamo.R, Bp_final, levels=levels_p, colors='k',
353 linewidths=0.8, alpha=0.4)
354 ax.set_xlabel('Latitude (degrees)')
355 ax.set_ylabel('Radius')
356 ax.set_title('Poloidal Field B_p (final)')
357 ax.axvline(x=0, color='k', linestyle='--', alpha=0.5)
358 plt.colorbar(cf, ax=ax, label='B_p')
359
360 # Toroidal field
361 ax = axes[1]
362 levels_t = np.linspace(-np.abs(Bt_final).max(), np.abs(Bt_final).max(), 21)
363 cf = ax.contourf(lat_2d, dynamo.R, Bt_final, levels=levels_t, cmap='RdBu_r', extend='both')
364 ax.contour(lat_2d, dynamo.R, Bt_final, levels=levels_t, colors='k',
365 linewidths=0.8, alpha=0.4)
366 ax.set_xlabel('Latitude (degrees)')
367 ax.set_ylabel('Radius')
368 ax.set_title('Toroidal Field B_t (final)')
369 ax.axvline(x=0, color='k', linestyle='--', alpha=0.5)
370 plt.colorbar(cf, ax=ax, label='B_t')
371
372 plt.tight_layout()
373 plt.savefig(f'/opt/projects/01_Personal/03_Study/examples/Numerical_Simulation/10_projects/{save_prefix}_fields.png',
374 dpi=150, bbox_inches='tight')
375 plt.close()
376 print(f"Figure saved: {save_prefix}_fields.png")
377
378
379def study_dynamo_regimes():
380 """
381 Study different dynamo regimes by varying dynamo numbers.
382 """
383 print("=" * 70)
384 print("Alpha-Omega Dynamo: Regime Study")
385 print("=" * 70)
386
387 fig, axes = plt.subplots(2, 2, figsize=(14, 10))
388
389 regimes = [
390 {'name': 'Subcritical', 'C_alpha': 0.5, 'C_omega': 2.0}, # D = 1.0
391 {'name': 'Critical', 'C_alpha': 1.0, 'C_omega': 3.0}, # D = 3.0
392 {'name': 'Supercritical', 'C_alpha': 1.5, 'C_omega': 5.0}, # D = 7.5
393 {'name': 'Strongly Supercritical', 'C_alpha': 2.0, 'C_omega': 8.0}, # D = 16
394 ]
395
396 for idx, regime in enumerate(regimes):
397 print(f"\n{regime['name']} Regime:")
398 print(f" C_α={regime['C_alpha']}, C_Ω={regime['C_omega']}")
399 print(f" Dynamo number D = {regime['C_alpha'] * regime['C_omega']}")
400
401 dynamo = AlphaOmegaDynamo(r_inner=0.3, r_outer=1.0, nr=40, ntheta=30)
402 dynamo.run_dynamo(t_final=80.0, dt=0.05,
403 C_alpha=regime['C_alpha'],
404 C_omega=regime['C_omega'],
405 eta=0.1)
406
407 # Extract time series
408 times = np.array(dynamo.time_history)
409 Bt_max = [np.max(np.abs(Bt)) for Bt in dynamo.Bt_history]
410
411 # Plot
412 ax = axes[idx // 2, idx % 2]
413 ax.plot(times, Bt_max, 'b-', linewidth=2)
414 ax.set_xlabel('Time')
415 ax.set_ylabel('|Bt|_max')
416 ax.set_title(f"{regime['name']}: D={regime['C_alpha']*regime['C_omega']:.1f}")
417 ax.grid(True, alpha=0.3)
418 ax.set_yscale('log')
419
420 # Analyze
421 if len(times) > 50:
422 final_avg = np.mean(Bt_max[-20:])
423 if final_avg < 1e-4:
424 result = "Decays (subcritical)"
425 elif np.std(Bt_max[-20:]) / np.mean(Bt_max[-20:]) > 0.3:
426 result = "Oscillatory (supercritical)"
427 else:
428 result = "Saturated"
429 print(f" Result: {result}")
430
431 plt.tight_layout()
432 plt.savefig('/opt/projects/01_Personal/03_Study/examples/Numerical_Simulation/10_projects/dynamo_regimes.png',
433 dpi=150, bbox_inches='tight')
434 plt.close()
435 print("\nFigure saved: dynamo_regimes.png")
436
437
438def main():
439 """Run complete alpha-omega dynamo study."""
440 print("=" * 70)
441 print("Alpha-Omega Dynamo in Spherical Shell")
442 print("=" * 70)
443
444 # Main simulation
445 print("\n[1] Running main dynamo simulation...")
446 dynamo = AlphaOmegaDynamo(r_inner=0.35, r_outer=1.0, nr=50, ntheta=40)
447 dynamo.run_dynamo(t_final=100.0, dt=0.05,
448 C_alpha=1.2, C_omega=6.0, eta=0.1)
449
450 # Visualize
451 visualize_dynamo(dynamo, save_prefix="dynamo_main")
452
453 # Study different regimes
454 print("\n[2] Studying different dynamo regimes...")
455 study_dynamo_regimes()
456
457 # Summary
458 print("\n" + "=" * 70)
459 print("Summary: Alpha-Omega Dynamo")
460 print("=" * 70)
461 print("""
462 Mean-Field Dynamo Theory:
463
464 The alpha-omega dynamo sustains magnetic fields through two mechanisms:
465
466 1. Ω-Effect (Differential Rotation):
467 - Shears poloidal field B_p into toroidal field B_t
468 - ∂B_t/∂t ∝ (∇Ω) · (∇B_p)
469 - Strong in regions with velocity gradients
470
471 2. α-Effect (Helical Turbulence):
472 - Twists toroidal field B_t back into poloidal field B_p
473 - ∂B_p/∂t ∝ α B_t
474 - Arises from cyclonic convection in rotating spheres
475
476 Critical Dynamo Number:
477 - D = C_α * C_Ω (dimensionless)
478 - D < D_crit: Field decays (subcritical)
479 - D ≈ D_crit: Marginal (critical dynamo number ~ 2-5)
480 - D > D_crit: Self-sustained oscillations (supercritical)
481
482 Butterfly Diagram:
483 - Shows latitudinal migration of magnetic activity
484 - In Sun: activity starts at ~30° latitude and migrates equatorward
485 - Period: ~11 years for solar cycle (22 years for full magnetic cycle)
486 - Our model shows similar equatorward migration pattern
487
488 Field Reversals:
489 - Poloidal field reverses periodically
490 - Earth's magnetic field reverses every ~200,000-300,000 years
491 - Mechanism: Nonlinear saturation and turbulent fluctuations
492
493 Applications:
494 - Solar dynamo: Sunspot cycle, solar flares
495 - Earth's dynamo: Geomagnetic field, pole reversals
496 - Stellar dynamos: Starspot activity, magnetic braking
497 - Galactic dynamo: Large-scale magnetic fields in galaxies
498
499 Limitations of This Model:
500 - Simplified geometry (spherical shell, not full 3D)
501 - Mean-field approximation (averages over turbulence)
502 - Linear α-effect (real α depends on field strength)
503 - No magnetic buoyancy or field instabilities
504 - Low spatial resolution
505
506 Advanced Topics:
507 - α²-dynamo: Both poloidal and toroidal generated by α
508 - Magnetic buoyancy: Rising flux tubes (sunspots)
509 - Tachocline: Interface dynamo at base of convection zone
510 - Grand minima: Maunder Minimum (1645-1715, reduced sunspots)
511 """)
512
513
514if __name__ == "__main__":
515 main()