"""Rikitake's reversal model, as derived by Janosi & Herein (2026): our re-run. Fact-Checked Science, October 2026, for the page "Rikitake's dynamo gets a physical derivation, and its reversals come from inside the core" (/rikitake-mean-field/). The paper: D. Janosi and M. Herein, "Grounding Rikitake-like conceptual dynamos in mean-field MHD: the derivation and geophysical scaling of the most basic reversal dynamo model", Geophysical Journal International 247(3), ggag390 (2026), https://doi.org/10.1093/gji/ggag390, CC BY 4.0. No code or data came with it; this script is our own, written from the equations it prints. Standard library only (Python 3.8 or later). It integrates the paper's canonical Rikitake model, eq. (41), x' = -x + y z, y' = -y - a x + x z, z' = c - x y, with the classical fourth-order Runge-Kutta method from x = 1, y = 0, z = a, the paper's initial conditions (caption of its Fig. 1). x and y are the amplitudes of two large-scale magnetic patterns, z that of one flow pattern; time is in units of t0 = 1/beta, the magnetic diffusion time, which the paper puts at 1.70e12 s, about 54,000 years (eq. 68). A polarity reversal is a change of sign of x. To ignore grazing of zero, a change is counted only when x passes from beyond +h to beyond -h or back, with h one thousandth of the size of x at the model's steady state, sqrt(c / (z* - a)) with z* = (a + sqrt(a^2 + 4)) / 2. A chron is the time between successive reversals. It prints five parts: 1 The paper's scalings recomputed from its inputs (eqs 50-80), and its conversions of model time into years in section 3.2, which come out 100 times too small. 2 Fig. 1 (a = 10, c = 1e4): the settled oscillation, over the 10 time units of a run of 10^3 steps of t0/100. 3 Fig. 2 (a = 8; c = 0.22, 0.31, 0.4, 0.6): the time of the first reversal at three step sizes, then reversals and chron lengths over 1,000 time units (10^5 steps of t0/100) and 10,000 (10^6 steps), in years at t0 = 53.95 kyr. 4 Where the model reverses: c from 0.001 to 1,000 at a = 8, and the harmonic realization's Earth-based range (a' = 0.082-8.2, c' = 10-1,000, eq. 99), each run classed by its behaviour after the first 100 time units as "irregular" (chron lengths vary, coefficient of variation 0.2 or more), "clockwork" (they don't) or "none" (no reversal). 5 With --large: the rest of the sweep, c from 1,780 to 562,000 (the paper's Earth-based range is 5,660-566,000, eq. 63), and long runs at three Earth-range settings. Run it as python3 rikitake_rerun.py parts 1-4, about a minute on 8 cores python3 rikitake_rerun.py --large all five parts, 5 minutes or more on 8 cores (its longest run cannot be split across cores) python3 rikitake_rerun.py --json out.json also writes the chart data to out.json Runs are spread over the computer's cores with the standard multiprocessing module; add --serial to use one core (about eight times slower). """ import json import math import multiprocessing as mp import statistics import sys YEAR = 365.25 * 86400.0 T0_KYR = 1.0 / (3.0 / 2.26e6 ** 2) / YEAR / 1e3 # t0 = L^2 / eta_T, about 53.95 kyr # ---------------------------------------------------------------- the model and its integration def integrate(a, c, dt, t_end, t_record=0.0, track=None): """RK4 from x=1, y=0, z=a; returns the times of the reversals of x up to t_end. track = (t_from, t_to) also returns the mean and the min/max of x over that interval.""" x, y, z = 1.0, 0.0, a zs = 0.5 * (a + math.sqrt(a * a + 4.0)) h = 1e-3 * math.sqrt(c / (zs - a)) h2, h6 = 0.5 * dt, dt / 6.0 n = int(round(t_end / dt)) state = 1 revs = [] s_sum, s_n, s_min, s_max = 0.0, 0, float("inf"), float("-inf") tr0, tr1 = track if track else (2.0 * t_end, 2.0 * t_end) for i in range(1, n + 1): k1x = -x + y * z; k1y = -y - a * x + x * z; k1z = c - x * y x2 = x + h2 * k1x; y2 = y + h2 * k1y; z2 = z + h2 * k1z k2x = -x2 + y2 * z2; k2y = -y2 - a * x2 + x2 * z2; k2z = c - x2 * y2 x3 = x + h2 * k2x; y3 = y + h2 * k2y; z3 = z + h2 * k2z k3x = -x3 + y3 * z3; k3y = -y3 - a * x3 + x3 * z3; k3z = c - x3 * y3 x4 = x + dt * k3x; y4 = y + dt * k3y; z4 = z + dt * k3z k4x = -x4 + y4 * z4; k4y = -y4 - a * x4 + x4 * z4; k4z = c - x4 * y4 x += h6 * (k1x + 2.0 * k2x + 2.0 * k3x + k4x) y += h6 * (k1y + 2.0 * k2y + 2.0 * k3y + k4y) z += h6 * (k1z + 2.0 * k2z + 2.0 * k3z + k4z) if x > h: if state < 0: revs.append(i * dt) state = 1 elif x < -h: if state > 0: revs.append(i * dt) state = -1 if track: t = i * dt if tr0 <= t <= tr1: s_sum += x; s_n += 1 if x < s_min: s_min = x if x > s_max: s_max = x if track: return revs, (s_sum / s_n, s_min, s_max) return revs def step_for(a, c): """The step we use: small enough for the fast motion at large c and large a.""" return min(0.005, 0.002 / math.sqrt(max(c, 1.0)), 0.02 / max(a, 1.0)) def classify(job): """One point of the sweep: (a, c, t_end) -> dict with the behaviour after t = 100.""" a, c, t_end = job dt = step_for(a, c) revs = integrate(a, c, dt, t_end) late = [t for t in revs if t > 100.0] out = {"a": a, "c": c, "dt": dt, "t_end": t_end, "n_reversals": len(revs), "first_reversal": revs[0] if revs else None, "n_after_100": len(late)} if len(late) >= 3: gaps = [t1 - t0 for t0, t1 in zip(late, late[1:])] mean = statistics.fmean(gaps) cv = statistics.pstdev(gaps) / mean out.update({"regime": "irregular" if cv >= 0.2 else "clockwork", "cv": cv, "median_chron": statistics.median(gaps), "longest_chron": max(gaps)}) elif late: out["regime"] = "rare" else: out["regime"] = "none" return out def fig2_run(c): """Fig. 2's settings, a = 8: reversal times over 10,000 time units at dt = 0.005.""" return c, integrate(8.0, c, 0.005, 10000.0) def first_reversal(job): c, dt = job revs = integrate(8.0, c, dt, 40.0) return c, dt, (revs[0] if revs else None) def chron_stats(revs, t_end): gaps = [t1 - t0 for t0, t1 in zip(revs, revs[1:])] if not gaps: return {"n_reversals": len(revs)} return {"n_reversals": len(revs), "shortest": min(gaps), "median": statistics.median(gaps), "longest": max(gaps), "first_interval": revs[0], "last_interval": t_end - revs[-1], "per_Myr": len(revs) / (t_end * T0_KYR / 1e3)} # ---------------------------------------------------------------- part 1: the paper's arithmetic def part1(): print("PART 1. The paper's scalings, recomputed from its inputs (section 3)") L, rho0, etaT, mu0 = 2.26e6, 1e4, 3.0, 4e-7 * math.pi # eqs 50-52 U = (1e-4, 1e-3) # eq. 56 alpha_e = (1e-6, 1e-4) # eq. 57 beta = etaT / L ** 2 # eq. 54 alpha = [ae / L for ae in alpha_e] # eq. 55 gamma = [u * u / L for u in U] # eq. 60 mu = 1.0 / L # eq. 62 a = [al / beta for al in alpha] c = [g * mu / beta ** 2 for g in gamma] delta = [beta / (4 * mu0 * rho0 * u) for u in U] # eq. 73 x0 = [beta / math.sqrt(mu * d) for d in delta] # eq. 76 rows = [ ("beta = eta_T / L^2 (eq. 54)", f"{beta:.3e} s^-1", "5.87e-13"), ("t0 = 1 / beta (eq. 68)", f"{1 / beta:.3e} s = {T0_KYR:.2f} kyr", "1.70e12 s, 54 kyr"), ("a = alpha / beta (eq. 59)", f"{a[0]:.3f} to {a[1]:.1f}", "0.749-74.9"), ("c = gamma mu / beta^2 (eq. 63)", f"{c[0]:.0f} to {c[1]:.0f}", "5.66e3-5.66e5"), ("Earth's c over the reversing c ~ 0.01-1 (eq. 64)", f"{c[0] / 1:,.0f} to {c[1] / 0.01:,.0f} times", "'an about five order of magnitude difference'"), ("delta (eq. 74)", f"{delta[1]:.3e} to {delta[0]:.3e}", "1.16e-8 to 1.16e-7"), ("x0 = beta / sqrt(mu delta) (eq. 76)", f"{x0[0] * 1e9:.0f} to {x0[1] * 1e9:.0f} nT", "2,500-8,000 nT"), ("z0 = beta / mu (eq. 77)", f"{beta / mu:.3e} m/s", "2.59e-7 (a slip; eta_T / L = 1.33e-6)"), ] for name, ours, paper in rows: print(f" {name:<50s} ours: {ours:<30s} paper: {paper}") print("\n The paper's conversions of model time into years (section 3.2), with dt = t0 / 100:") step_yr = T0_KYR * 1e3 / 100 for label, steps, paper in (("Fig. 1's run, 10^3 dt", 1e3, "5.4 kyr"), ("Fig. 2's shorter runs, 10^5 dt", 1e5, "540 kyr"), ("Fig. 2's longer run, 10^6 dt", 1e6, "5.4 Myr")): yrs = steps * step_yr txt = f"{yrs / 1e3:,.0f} kyr" if yrs < 1e6 else f"{yrs / 1e6:,.0f} Myr" print(f" {label:<34s} = {steps * 0.01:>7,.0f} t0 = {txt:<10s} the paper says {paper}") print(" The longest chrons it prints, back in model time at t0 = 54 kyr:") for cv, kyr in ((0.22, 54), (0.31, 1620), (0.4, 54), (0.6, 16.2)): print(f" c = {cv}: {kyr:>6} kyr = {kyr / T0_KYR:5.2f} time units") return {"t0_kyr": T0_KYR, "c_earth": c, "a_earth": a, "x0_nT": [v * 1e9 for v in x0]} # ---------------------------------------------------------------- part 2: Fig. 1 def part2(): print("\nPART 2. Fig. 1 (a = 10, c = 1e4): the settled oscillation") res = {} for t_end in (0.1, 10.0): dt = step_for(10.0, 1e4) / 2 revs, (m, lo, hi) = integrate(10.0, 1e4, dt, t_end, track=(t_end / 2, t_end)) res[str(t_end)] = {"mean": m, "half_range": (hi - lo) / 2, "reversals": len(revs)} print(f" run of {t_end:>4} units: second half mean of x {m:6.1f}, half-range {(hi - lo) / 2:6.1f}," f" reversals {len(revs)}") print(" The paper's caption: 'a quasi-periodic oscillation with an amplitude of about 150 units around" " the value 450' after 'around 3 units'.") return res # ---------------------------------------------------------------- part 3: Fig. 2 def part3(pool_map): print("\nPART 3. Fig. 2 (a = 8): first reversals, then chrons over 1,000 and 10,000 time units") cs = (0.22, 0.31, 0.4, 0.6) firsts = pool_map(first_reversal, [(c, dt) for c in cs for dt in (0.01, 0.005, 0.0025)]) for c in cs: vals = [f"{fr:.3f} (dt {dt})" for cc, dt, fr in firsts if cc == c] print(f" c = {c}: first reversal at t = " + ", ".join(vals)) runs = dict(pool_map(fig2_run, cs)) out = {"t0_kyr": T0_KYR, "runs": []} for c in cs: revs = runs[c] r1 = [t for t in revs if t <= 1000.0] s1, s10 = chron_stats(r1, 1000.0), chron_stats(revs, 10000.0) out["runs"].append({"c": c, "reversals_1000": r1, "stats_1000": s1, "stats_10000": s10, "first_reversal": revs[0]}) for label, s, span in (("1,000 units (54 Myr)", s1, 1000.0), ("10,000 units (540 Myr)", s10, 10000.0)): if "longest" not in s: print(f" c = {c}, {label}: {s['n_reversals']} reversals") continue print(f" c = {c}, {label}: {s['n_reversals']:4d} reversals ({s['per_Myr']:.2f} per Myr);" f" chrons: shortest {s['shortest']:6.2f} units = {s['shortest'] * T0_KYR:6.0f} kyr," f" median {s['median']:6.2f} = {s['median'] * T0_KYR:6.0f} kyr," f" longest {s['longest']:8.2f} = {s['longest'] * T0_KYR / 1e3:7.2f} Myr") out["first_reversal_checks"] = [{"c": c, "dt": dt, "t": fr} for c, dt, fr in firsts] return out # ---------------------------------------------------------------- parts 4 and 5: the sweep SWEEP_SMALL = [round(10 ** (k / 4.0), 6) for k in range(-12, 13)] # 0.001 ... 1000 SWEEP_LARGE = [round(10 ** (k / 4.0), 1) for k in range(13, 24)] # 1778 ... 562341 HARMONIC = [(a, c) for a in (0.082, 0.82, 8.2) for c in (10.0, 100.0, 1000.0)] EARTH_LONG = [(10.0, 1e4), (7.49, 5660.0), (0.749, 5660.0)] def t_end_for(c): return 2100.0 if c <= 10 else 1100.0 def show(r, t0_kyr=T0_KYR): first = f"{r['first_reversal']:.2f}" if r["first_reversal"] is not None else "-" if r["regime"] in ("irregular", "clockwork"): med = r["median_chron"] when = f"{med * t0_kyr:,.0f} kyr" if med * t0_kyr >= 10 else f"{med * t0_kyr * 1e3:,.0f} years" return (f" a = {r['a']:<6g} c = {r['c']:<10g} {r['regime']:<9s} first at {first:>8s};" f" median chron {med:9.4f} units = {when:>11s}; CV {r['cv']:.3f}; {r['n_after_100']} after t = 100") return f" a = {r['a']:<6g} c = {r['c']:<10g} {r['regime']:<9s} first at {first:>8s}; {r['n_after_100']} after t = 100" def part4(pool_map, large): print("\nPART 4. Where the model reverses, at a = 8 (the value of Fig. 2), and in the harmonic" " realization's Earth-based range") jobs = [(8.0, c, t_end_for(c)) for c in SWEEP_SMALL] + [(a, c, t_end_for(c)) for a, c in HARMONIC] if large: jobs += [(8.0, c, t_end_for(c)) for c in SWEEP_LARGE] + [(a, c, 1100.0) for a, c in EARTH_LONG] res = pool_map(classify, jobs) sweep = sorted([r for r in res if r["a"] == 8.0 and (r["c"] in SWEEP_SMALL or r["c"] in SWEEP_LARGE)], key=lambda r: r["c"]) harm = [r for r in res if (r["a"], r["c"]) in HARMONIC] earth = [r for r in res if (r["a"], r["c"]) in EARTH_LONG] print(" The sweep at a = 8:") for r in sweep: print(show(r)) t0h_kyr = T0_KYR / (4 * math.pi ** 2) # t0' = 1/beta' = L^2 / (4 pi^2 eta_T), eqs 93 and 100 print(" The harmonic realization's Earth-based range, a' = 0.082-8.2 and c' = 10-1,000 (eq. 99);" f" its time unit is t0' = 1/beta' = {t0h_kyr:.2f} kyr (eq. 100):") for r in harm: print(show(r, t0h_kyr)) if large: print("\nPART 5. Long runs at the paper's Earth-range settings (1,100 time units, 59 Myr)") for r in earth: print(show(r)) return {"sweep_a8": sweep, "harmonic": harm, "earth_long": earth} def main(): args = sys.argv[1:] large = "--large" in args serial = "--serial" in args out_path = args[args.index("--json") + 1] if "--json" in args else None result = {"part1": part1(), "part2": part2()} if serial: pool_map = lambda f, jobs: [f(j) for j in jobs] result["part3"] = part3(pool_map) result["part4"] = part4(pool_map, large) else: with mp.Pool() as pool: pool_map = lambda f, jobs: pool.map(f, jobs, chunksize=1) result["part3"] = part3(pool_map) result["part4"] = part4(pool_map, large) if out_path: with open(out_path, "w", encoding="utf-8") as fh: json.dump(result, fh, indent=1) print(f"\nWrote {out_path}") if __name__ == "__main__": main()