Summary

Key papers (10)

1. Is there a chronometer hidden deep in the Sun?. R. H. Dicke; Nature 276, 676-680, 1978.
   Poses the clock-versus-random-walk question and defines Dicke's ratio with its two limiting expectations; our study is a direct, modern re-execution of this test on the recalibrated record.
2. Is the solar cycle timed by a clock?. P. Hoyng; Solar Physics 169, 253-264, 1996.
   Shows that alpha-fluctuations in mean-field dynamo theory reproduce the apparent phase memory, so no extra timing device is required either observationally or theoretically.
3. A Model of a Tidally Synchronized Solar Dynamo. F. Stefani, A. Giesecke, T. Weier; Solar Physics 294, 60, 2019.
   The principal modern claim that the cycle is clocked, at the 11.07-yr Venus-Earth-Jupiter alignment period, argued from Dicke's ratio on a 1000-yr epoch series.
4. Tidally Synchronized Solar Dynamo: A Rebuttal. H.-C. Nataf; Solar Physics 297, 107, 2022.
   Shows the epochs underpinning the clock claim are clocked by construction because Schove assumed a continued 11-yr cycle; also states the Dicke-ratio formulae explicitly.
5. No Evidence for Absence of Solar Dynamo Synchronization. F. Stefani, J. Beer, T. Weier; Solar Physics 298, 83, 2023.
   The counter-reply, defending a stable 11.07-yr period from cosmogenic-radionuclide reconstructions; documents that the controversy is unresolved and still rests on proxy epochs.
6. Dynamo models of the solar cycle. P. Charbonneau; Living Reviews in Solar Physics 17, 4, 2020.
   Reference account of stochastic forcing of the alpha-effect and Babcock-Leighton source, the theoretical setting in which a diffusing cycle phase is expected.
7. Solar activity: periodicities beyond 11 years are consistent with random forcing. R. H. Cameron, M. Schuessler; Astronomy & Astrophysics 625, A28, 2019.
   Demonstrates that apparent long-period regularities in solar activity are realization noise of a noisy limit cycle, the null against which any claim of extra timekeeping must be judged.
8. A history of solar activity over millennia. I. G. Usoskin; Living Reviews in Solar Physics 20, 2, 2023.
   Sets out the reliability limits of Schove-type epoch reconstructions, motivating a test restricted to directly measured sunspot numbers.
9. Revisiting the Sunspot Number. A 400-Year Perspective on the Solar Cycle. F. Clette, L. Svalgaard, J. M. Vaquero, E. W. Cliver; Space Science Reviews 186, 35-103, 2014.
   Documents the recalibration behind SILSO version 2.0, whose values are not comparable with the version 1.0 numbers used in all pre-2015 analyses of this question.
10. The Solar Cycle. D. H. Hathaway; Living Reviews in Solar Physics 12, 4, 2015.
    Standard compilation of cycle-length statistics, the Waldmeier effect and the length-amplitude anticorrelation, providing the descriptive baseline for any phase model.

1. Idea 1Selected

   No chronometer in the Sun: the Schwabe cycle phase diffuses as a random walk

   Hypothesis. The epochs of sunspot-cycle minima are NOT anchored to a rigid clock. Writing the n-th minimum as t_n = t_0 + nP + xi_n + W_n, with xi_n an independent rigid-clock jitter of rms sigma_c and W_n = sum_{k<=n} eps_k a random walk of the phase with per-cycle step rms sigma_w, the hypothesis is sigma_w > 0 with the clock limit sigma_w = 0 rejected. Two parameter-free corollaries follow and are tested as well: the lag-1 autocorrelation of successive cycle lengths, which a rigid clock forces to exactly -1/2, must instead be consistent with 0; and Dicke's ratio Di(25) must lie near its random-walk expectation (n+2)(n^2-1)/(3(5n^2+6n-3)) rather than its clock expectation (n^2-1)/(2(n^2+2n+3)). The hypothesis is refuted if sigma_w is consistent with zero, or if the likelihood-ratio test fails to reject sigma_w = 0 at the 5 per cent level, or if the lag-1 autocorrelation is consistent with -1/2.

   Dataset. SILSO monthly mean total sunspot number, version 2.0 sidc_sunspot_monthly

   Test plan, novelty and risk

   Test plan. VARIABLES: decimal_year, sunspot_number, sd and n_observations from the monthly file. EPOCHS: smooth sunspot_number with the SIDC 13-month kernel (unit weights with half weight on the two end months); take every month at which the smoothed series attains its minimum over a +-48-month window, merging months within 12 of one another into one epoch and using their mean date, which resolves the flat 1810 minimum. This yields the cycle-1 to cycle-25 minima, i.e. 25 epochs and 24 cycle lengths. Cycle maxima are located the same way as an independent marker. MEASUREMENT UNCERTAINTY: the standard error of each monthly mean is sd/sqrt(n_observations), available from 1818; draw 600 Monte Carlo replicates of the monthly series with independent Gaussian perturbations of that size, re-smooth and re-locate every epoch, and take the replicate standard deviation as the epoch measurement uncertainty s_n, adding in quadrature the variance (1/12)^2/12 of the one-month quantisation of an epoch. SAMPLE FILTERS: primary sample all 25 minima; a secondary sample restricted to the 19 minima from 1823 onward, for which every contributing month carries a published dispersion. ESTIMATOR: Gaussian maximum likelihood for the covariance Cov(t_n,t_m) = sigma_c^2 delta_nm + sigma_w^2 min(n,m) + s_n^2 delta_nm with the mean t_0 + nP profiled out by generalised least squares; report P and its standard error, sigma_w with a profile-likelihood 68 per cent interval and sigma_c with a 95 per cent upper limit. TESTS: likelihood ratio for sigma_w = 0 (a boundary case, so the asymptotic p-value uses the 50:50 mixture of chi^2_0 and chi^2_1) and, because N is small, an exact p-value from 2000 parametric-bootstrap replicates simulated under the maximum-likelihood clock; likelihood ratio for sigma_c = 0; the lag-1 autocorrelation of the 24 cycle lengths against 20000 simulations of each model; Dicke's ratio against the same simulations and against both analytic expectations; and the phase structure function S(k) = <(t_{n+k}-t_n-kP)^2> for k = 1..12, which must grow linearly under a random walk and saturate under a clock. DERIVED QUANTITIES: phase-diffusion coefficient sigma_w^2/P^3 in squared cycles per year, the coherence horizon (P/sigma_w)^2 cycles, the residuals of the epochs against a rigid 11.07-yr planetary-tidal ephemeris, and a falsifiable forecast of the cycle-26 minimum with its uncertainty. ROBUSTNESS ACTUALLY RUN: maxima instead of minima; the 1823-onward subsample; a 24-month smoothing kernel; refitting with the measurement uncertainties set to zero to size their contribution; and a quadratic mean term to test for a secular drift of the period. REFUTATION: sigma_w consistent with zero, or the clock not rejected, or lag-1 autocorrelation consistent with -1/2.

   Novelty. The strongest published claim that the cycle is clocked (Stefani, Giesecke and Weier 2019) needs a millennium of epochs to reach significance, and those pre-telescopic epochs come from Schove's reconstruction, which Nataf (2022) showed is clocked by construction because Schove assumed a continued 11-yr cycle. Nobody has settled the question on directly measured sunspot numbers alone, and nobody in the exchange propagates the published per-month measurement uncertainty into the epochs. Fitting a model that nests both hypotheses, rather than choosing between two point hypotheses, and delivering a measured phase-diffusion coefficient with an uncertainty, are both new.

   Expected outcome. A per-cycle phase step of order 1 yr, significantly non-zero, with the rigid-clock component consistent with zero and an upper limit of a few tenths of a year; the clock model rejected at better than the 1 per cent level; lag-1 autocorrelation near 0 rather than -1/2; Dicke's ratio near its random-walk expectation; and a coherence horizon of order 100 cycles, i.e. about a millennium, after which the cycle phase is no longer predictable.

   Main risk. Only about two dozen cycles exist, so the power is limited and the confidence intervals on sigma_w will be wide; and a slow deterministic drift of the period, of Gleissberg type, would also inflate the long-lag phase dispersion and cannot be cleanly separated from genuine diffusion with 24 cycles. Both must be stated explicitly rather than glossed over.

2. Idea 2

   Sub-proportional scaling of the seasonal CO2 amplitude with the atmospheric CO2 burden

   Hypothesis. The peak-to-trough amplitude A(t) of the seasonal cycle of atmospheric CO2 at Mauna Loa scales with the trend mole fraction C(t) as A = A_0 (C/C_0)^gamma with gamma = 1, the value expected if the seasonal biospheric exchange is a fixed fraction of the ambient burden. The hypothesis is refuted if the maximum-likelihood gamma differs from 1 by more than twice its standard error.

   Dataset. Mauna Loa monthly mean atmospheric CO2 mauna_loa_co2_monthly

   Test plan, novelty and risk

   Test plan. VARIABLES: decimal_date, average and unc. SAMPLE FILTER: months from May 1974 onward, for which NOAA publishes unc, and excluding months with missing unc, which flag interpolated values; report the number of months actually used. MODEL: C(t) = T(t) + A_0 (T(t)/C_0)^gamma S(t) with T(t) a cubic polynomial trend and S(t) = sum_{k=1..4} [a_k cos(2 pi k t) + b_k sin(2 pi k t)] normalised to unit peak-to-trough amplitude. ESTIMATOR: weighted nonlinear least squares with weights 1/unc^2, gamma and its standard error from the Jacobian; because the residual scatter will exceed the stated measurement uncertainty, report both the measurement-limited error and an error inflated by the square root of the reduced chi-square, and a moving-block bootstrap error that respects the residual autocorrelation. TESTS: Wald test of gamma = 1; likelihood-ratio test of the nested constant-amplitude model gamma = 0. ROBUSTNESS: two and six harmonics; quadratic and quartic trends; splitting the record in half.

   Novelty. The growth of the Mauna Loa seasonal amplitude is documented, but it is almost always reported as a linear trend in years. Casting it as a scaling exponent in the CO2 burden itself gives a dimensionless, physically interpretable number with a null of exactly 1, and makes the measurement-uncertainty budget explicit.

   Expected outcome. A gamma significantly different from 1, with the reduced chi-square far above 1, showing that interannual geophysical variability rather than measurement error dominates the uncertainty budget.

   Main risk. Trend and amplitude are strongly correlated in the fit, so gamma may be poorly identified; the trend polynomial order can absorb part of the amplitude change; and the residuals are strongly autocorrelated, so naive standard errors would be far too small.

3. Idea 3

   Which elemental property controls the critical temperature? A collinearity-resolved ranking

   Hypothesis. Of the eight elemental properties summarised in the SuperCon feature set, the stoichiometrically weighted statistics of valence and of thermal conductivity carry more independent predictive information about the measured critical temperature than those of atomic mass, so that grouping the 80 predictors into eight physical blocks and dropping each block in turn reduces the cross-validated coefficient of determination most for valence and thermal conductivity. The hypothesis is refuted if the atomic-mass block, which an isotope-effect argument would privilege, causes a drop at least as large as either of them.

   Dataset. Critical temperatures of 21,263 superconductors superconductor_tc

   Test plan, novelty and risk

   Test plan. VARIABLES: all 81 predictors and critical_temp. SAMPLE FILTER: rows with critical_temp above 0 K; report the count. ESTIMATOR: ridge regression and gradient-boosted trees on standardised predictors, both with the penalty or depth chosen by inner cross-validation, under a fixed-seed 10-fold outer cross-validation; the statistic is the mean drop in out-of-fold R^2 when one physical block of ten summary statistics is removed, with its standard error across the ten folds. Also report block-permutation importances. TESTS: paired t-tests across folds between the valence, thermal-conductivity and atomic-mass blocks.

   Novelty. The original benchmark paper reports single-feature importances from one model; grouping the extremely collinear predictors into physically meaningful blocks and measuring block-level drop-out importance with fold-level errors gives an interpretable answer to which elemental property actually matters.

   Expected outcome. Valence and thermal conductivity dominate, with atomic mass contributing little once the other blocks are present, consistent with the cuprate-dominated tail of the sample.

   Main risk. Blocks are mutually collinear, so drop-out importance understates a block whose information is duplicated elsewhere; and the sample is unevenly distributed over chemistry, so the ranking may describe SuperCon's sampling rather than physics.

Ranking rationale

Idea 1 is first because it settles a specific, long-standing and currently contested physical question with a single well-posed likelihood test, the null is parameter-free rather than fitted, the dataset is exactly the right one and is the only directly measured record in existence, and the measurement uncertainty the venue expects can be propagated from published per-month dispersions rather than assumed. It also yields a new measured constant, the phase-diffusion coefficient, and a falsifiable forecast. Idea 2 is a clean uncertainty-weighted fit with a sharp null but the underlying amplitude growth is already well documented, so the contribution is a reparameterisation rather than a new answer. Idea 3 is the weakest for this venue: it is a predictive-modelling exercise whose conclusion is about a database's sampling as much as about physics, and it has no propagated measurement uncertainty because the SuperCon critical temperatures carry none.

Data provenance

Result (RESULT_JSON)

{
  "n_months_total": 3332,
  "n_months_missing_ssn": 0,
  "n_months_with_stated_error": 2503,
  "n_months_provisional": 5,
  "record_start": 1749.042,
  "record_end": 2026.623,
  "median_monthly_se": 1.1502173707608487,
  "median_monthly_se_1818_1849": 2.193931022920578,
  "median_monthly_se_2000_2026": 0.23012177511630372,
  "n_minima": 25,
  "n_maxima": 24,
  "first_minimum": 1755.123,
  "last_minimum": 2019.958,
  "epoch_sigma_median": 0.024056261216234408,
  "epoch_sigma_max": 0.24779782983082185,
  "epoch_sigma_max_year": 1902.042,
  "n_epochs_unmoved_by_noise": 14,
  "epoch_quantisation_sigma": 0.024056261216234408,
  "n_cycles": 24,
  "P": 11.034786928655784,
  "P_se": 0.23570072202057415,
  "mean_cycle_length": 11.034791666666669,
  "sd_cycle_length": 1.1793288716725157,
  "shortest_cycle": 9,
  "longest_cycle": 13.579999999999927,
  "sigma_w": 1.1546721277171903,
  "sigma_w_lo68": 1.0060378369054217,
  "sigma_w_hi68": 1.3435620081549022,
  "sigma_c": 0,
  "sigma_c_ul95": 0.15637456445529277,
  "sigma_c_under_clock": 1.6652504040922345,
  "LR_clock": 27.04321810856809,
  "p_clock_asymptotic": 9.947851308793453e-8,
  "p_clock_bootstrap": 0.0004997501249375312,
  "n_bootstrap_LR": 2000,
  "significance_sigma_clock": 5.200309424310066,
  "LR_randomwalk": 6.323830348264892e-12,
  "p_randomwalk": 0.4999989967707843,
  "r1": 0.10270290807450956,
  "r1_lo68": -0.09898407872441536,
  "r1_hi68": 0.3068498108917252,
  "r1_clock_expectation": -0.5,
  "r1_clock_sim_mean": -0.47903529974407455,
  "r1_clock_sim_sd": 0.14687010122711505,
  "p_r1_vs_clock": 0.0008499575021248938,
  "r1_rw_sim_mean": -0.0477155151001121,
  "r1_rw_sim_sd": 0.199161197338526,
  "p_r1_vs_randomwalk": 0.4617269136543173,
  "dicke_ratio": 2.171341196744508,
  "dicke_n": 25,
  "dicke_expect_clock": 0.46017699115044247,
  "dicke_expect_randomwalk": 1.7163814180929096,
  "dicke_clock_sim_mean": 0.4787540407758899,
  "dicke_clock_sim_sd": 0.09851936235252975,
  "dicke_rw_sim_mean": 1.6581351936870905,
  "dicke_rw_sim_sd": 0.8930516341714142,
  "p_dicke_vs_clock": 0.00004999750012499375,
  "p_dicke_vs_randomwalk": 0.4217289135543223,
  "dicke_ratio_maxima": 1.0667854187107084,
  "sf_kmax": 12,
  "sf_slope": 0.9136315947685262,
  "sf_slope_sim_sd": 0.6802775832435396,
  "sf_intercept": 0.9792681452962064,
  "sf_clock_sim_slope_mean": 0.08617729533148583,
  "sf_clock_sim_slope_sd": 0.21775059602422683,
  "sf_rw_sim_slope_mean": 0.643303581307581,
  "p_sf_vs_clock": 0.00899775056235941,
  "n_bootstrap_LR_exceeding": 0,
  "endpoint_period": 11.034791666666669,
  "sigma_w_no_meas_err": 1.1540893536228176,
  "meas_var_fraction_of_sigma_w2": 0.0036903534360983262,
  "maxima_n_cycles": 23,
  "maxima_P": 10.99274779332274,
  "maxima_P_se": 0.4081718361529133,
  "maxima_sigma_w": 1.9575105069781504,
  "maxima_sigma_c": 0,
  "maxima_LR_clock": 10.882845224876547,
  "maxima_p_clock": 0.00048529515669544986,
  "post1818_n_cycles": 18,
  "post1818_P": 10.921649235754396,
  "post1818_P_se": 0.2016126712997723,
  "post1818_sigma_w": 0.8551955469350631,
  "post1818_sigma_c": 0,
  "post1818_LR_clock": 17.132207754543813,
  "post1818_p_clock": 0.000017433006711023756,
  "smooth24_P": 10.992974885856416,
  "smooth24_sigma_w": 1.092065517871609,
  "smooth24_sigma_c": 0,
  "smooth24_LR_clock": 28.30382383366307,
  "smooth24_p_clock": 5.1845477332818345e-8,
  "LR_period_drift": 0.25084449733397207,
  "p_period_drift": 0.6164810676983199,
  "clock_resid_rms": 1.6691450564260175,
  "clock_resid_max_abs": 3.982029743588328,
  "clock_resid_max_year": 1784.708,
  "tidal_period": 11.07,
  "tidal_resid_rms": 1.669146692654662,
  "tidal_resid_max_abs": 3.984946666666474,
  "P_minus_tidal_sigma": -0.14939738428608948,
  "pred_cycle26_minimum": 2030.992786928656,
  "pred_cycle26_minimum_se": 1.1787286611394785,
  "phase_drift_per_cycle_deg": 37.67013977394716,
  "coherence_cycles": 91.3293860666246,
  "coherence_years": 1007.8003155701468,
  "phase_diffusion_cycles2_per_year": 0.0009922600584166973,
  "hypothesis_supported": true,
  "sf_obs_k1": 1.3328658964120204,
  "sf_rw_k1": 1.2886595224218345,
  "sf_clock_k1": 5.601982407648848,
  "sf_obs_k2": 3.0629651977656978,
  "sf_rw_k2": 2.454306989567359,
  "sf_clock_k2": 5.63778368598131,
  "sf_obs_k4": 4.468028895502645,
  "sf_rw_k4": 4.443056046743987,
  "sf_clock_k4": 5.660468372172896,
  "sf_obs_k8": 8.805153241829899,
  "sf_rw_k8": 7.118236982047629,
  "sf_clock_k8": 6.012154311954938,
  "sf_obs_k12": 10.190768438034096,
  "sf_rw_k12": 8.046949109396529,
  "sf_clock_k12": 6.573400815955307
}

"""Clock or random walk? Phase memory of the solar activity cycle.

Epochs of sunspot-cycle minima and maxima are located in the 13-month smoothed
SILSO v2.0 monthly total sunspot number and their measurement uncertainty is
propagated by Monte Carlo from the per-month standard error sd/sqrt(n_obs).
A phase model with rigid-clock jitter (sigma_c) plus a random walk (sigma_w) is
fitted by maximum likelihood; both limits are tested by likelihood ratio.
"""
import json
import numpy as np
from scipy import optimize, stats
from dataset import load

SEED, B_MC, B_SIM, B_LR = 20260908, 600, 20000, 2000
QV = (1.0 / 12.0) ** 2 / 12.0     # variance of one-month rounding of an epoch, yr^2
P_TIDAL = 11.07                   # period of the proposed planetary-tidal clock, yr

df = load().sort_values("decimal_year").reset_index(drop=True)
t = df["decimal_year"].to_numpy(float)
x = df["sunspot_number"].to_numpy(float)
sd, nob, yr = (df[c].to_numpy(float) for c in ("sd", "n_observations", "year"))
has_err = np.isfinite(sd) & np.isfinite(nob) & (nob > 0)
se = np.zeros_like(x)
se[has_err] = sd[has_err] / np.sqrt(nob[has_err])

W13 = np.ones(13); W13[0] = W13[-1] = 0.5; W13 /= W13.sum()
W24 = np.ones(25); W24[0] = W24[-1] = 0.5; W24 /= W24.sum()

def smooth(y, w):
    h = len(w) // 2
    s = np.convolve(y, w, mode="same"); s[:h] = np.nan; s[-h:] = np.nan
    return s

def epochs_of(s, half, kind):
    hits = [i for i in range(half, len(s) - half)
            if not np.isnan(s[i - half:i + half + 1]).any()
            and s[i] == (s[i - half:i + half + 1].min() if kind == "min"
                         else s[i - half:i + half + 1].max())]
    groups, cur = [], [hits[0]]
    for i in hits[1:]:
        if i - cur[-1] <= 12:
            cur.append(i)
        else:
            groups.append(cur); cur = [i]
    groups.append(cur)
    return np.array([float(np.mean(t[g])) for g in groups])

def relocate(s, centre, kind, halfwin=30):
    i0 = int(np.argmin(np.abs(t - centre)))
    lo, hi = max(6, i0 - halfwin), min(len(t) - 6, i0 + halfwin + 1)
    seg = s[lo:hi]
    v = np.nanmin(seg) if kind == "min" else np.nanmax(seg)
    return float(np.mean(t[np.where(seg == v)[0] + lo]))

def mc_sigma(nominal, kind, w, seed):
    rng = np.random.default_rng(seed)
    d = np.empty((B_MC, len(nominal)))
    for b in range(B_MC):
        s = smooth(x + rng.normal(0.0, se), w)
        d[b] = [relocate(s, c, kind) for c in nominal]
    return d.std(axis=0, ddof=1)

s13, s24 = smooth(x, W13), smooth(x, W24)
ep_min, ep_max = epochs_of(s13, 48, "min"), epochs_of(s13, 48, "max")
ep_min24 = np.array([relocate(s24, c, "min") for c in ep_min])
sig_min = mc_sigma(ep_min, "min", W13, SEED)
sig_max = mc_sigma(ep_max, "max", W13, SEED + 1)
sig_m24 = mc_sigma(ep_min24, "min", W24, SEED + 2)

# ------------------------------------------------------------- phase model
def make_nll(tt, vmeas, deg=1):
    N = len(tt); n = np.arange(N, dtype=float)
    X = np.column_stack([n ** k for k in range(deg + 1)])
    M = np.minimum.outer(np.arange(N), np.arange(N)).astype(float)
    def nll(vc, vw, want=False):
        C = max(vc, 0) * np.eye(N) + max(vw, 0) * M + np.diag(vmeas)
        try:
            L = np.linalg.cholesky(C)
        except np.linalg.LinAlgError:
            return (1e12, None) if want else 1e12
        z = np.linalg.solve(L, X); y = np.linalg.solve(L, tt)
        b, *_ = np.linalg.lstsq(z, y, rcond=None); r = y - z @ b
        v = float(np.sum(np.log(np.diag(L))) + 0.5 * float(r @ r))
        return (v, (b, np.linalg.inv(z.T @ z))) if want else v
    return nll, N

def mle(tt, vmeas, deg=1):
    nll, N = make_nll(tt, vmeas, deg)
    best = min((optimize.minimize(lambda p: nll(p[0], p[1]), p0, method="L-BFGS-B",
                                  bounds=[(0, 60), (0, 60)])
                for p0 in [(.5, .01), (.01, .5), (.3, .3), (1., 1e-9), (1e-9, 1.)]),
               key=lambda r: r.fun)
    vc, vw = float(best.x[0]), float(best.x[1])
    _, (beta, cb) = nll(vc, vw, want=True)
    fixw = lambda v: optimize.minimize_scalar(lambda u: nll(u, v), bounds=(0, 60),
                                              method="bounded")            # vw fixed
    fixc = lambda v: optimize.minimize_scalar(lambda u: nll(v, u), bounds=(0, 60),
                                              method="bounded")            # vc fixed
    return dict(N=N, vc=vc, vw=vw, nll=float(best.fun), beta=beta, cb=cb,
                fixw=fixw, fixc=fixc, nll_fn=nll)

def prof_ci(prof, v_mle, nll_min, delta, hi=60.0):
    g = lambda v: 2.0 * (prof(v).fun - nll_min) - delta
    out = []
    for a, b in ((0.0, v_mle), (v_mle, hi)):
        if b <= a or g(a) * g(b) > 0:
            out.append(float(np.sqrt(max(a if g(a) < 0 else b, 0.0))))
        else:
            out.append(float(np.sqrt(max(optimize.brentq(g, a, b, xtol=1e-12), 0.0))))
    return out

def analyse(tt, sig):
    vmeas = sig ** 2 + QV
    f = mle(tt, vmeas)
    r_clock, r_rw = f["fixw"](0.0), f["fixc"](0.0)              # sigma_w=0 ; sigma_c=0
    LRc, LRr = 2 * (r_clock.fun - f["nll"]), 2 * (r_rw.fun - f["nll"])
    sw_lo, sw_hi = prof_ci(f["fixw"], f["vw"], f["nll"], 1.0)
    _, sc_ul = prof_ci(f["fixc"], f["vc"], f["nll"], 3.8415)
    L = np.diff(tt)
    return dict(N=f["N"], ncyc=int(len(L)), P=float(f["beta"][1]),
                P_se=float(np.sqrt(f["cb"][1, 1])), sw=float(np.sqrt(f["vw"])),
                sw_lo=sw_lo, sw_hi=sw_hi, sc=float(np.sqrt(f["vc"])), sc_ul95=sc_ul,
                LR_clock=float(LRc), p_clock=0.5 * float(stats.chi2.sf(max(LRc, 0), 1)),
                LR_rw=float(LRr), p_rw=0.5 * float(stats.chi2.sf(max(LRr, 0), 1)),
                vc0=float(max(r_clock.x, 0.0)), meanL=float(L.mean()),
                sdL=float(L.std(ddof=1)), r1=float(np.corrcoef(L[:-1], L[1:])[0, 1]),
                nll=f["nll"], vmeas=vmeas)

main = analyse(ep_min, sig_min)
no_err = analyse(ep_min, np.zeros_like(sig_min))
maxima = analyse(ep_max, sig_max)
post = analyse(ep_min[6:], sig_min[6:])
sm24 = analyse(ep_min24, sig_m24)
LR_drift = 2.0 * (main["nll"] - mle(ep_min, main["vmeas"], deg=2)["nll"])

# ------------------------------------------------ simulation of both models
rng = np.random.default_rng(SEED + 7)
NC = main["ncyc"]; n_ax = np.arange(NC + 1, dtype=float)
SC0 = float(np.sqrt(main["vc0"]))          # ML clock amplitude under sigma_w = 0
def sim(sc, sw, nb):
    e = rng.normal(0.0, sw, (nb, NC + 1)) if sw > 0 else np.zeros((nb, NC + 1))
    e[:, 0] = 0.0
    return (main["P"] * n_ax + np.cumsum(e, axis=1) + rng.normal(0.0, max(sc, 1e-12), (nb, NC + 1))
            + rng.normal(0.0, np.sqrt(main["vmeas"]), (nb, NC + 1)))
def r1_of(T):
    L = np.diff(T, axis=1); a, b = L[:, :-1], L[:, 1:]
    az, bz = a - a.mean(1, keepdims=True), b - b.mean(1, keepdims=True)
    return (az * bz).sum(1) / np.sqrt((az ** 2).sum(1) * (bz ** 2).sum(1))
def two_sided(sample, obs):
    n = len(sample)
    return float((2 * min((sample <= obs).sum(), (sample >= obs).sum()) + 1) / (n + 1))
r1_ck, r1_rw = r1_of(sim(SC0, 0.0, B_SIM)), r1_of(sim(main["sc"], main["sw"], B_SIM))
p_r1_ck, p_r1_rw = two_sided(r1_ck, main["r1"]), two_sided(r1_rw, main["r1"])
r1_lo, r1_hi = [float(v) for v in np.percentile(r1_rw - r1_rw.mean() + main["r1"], [16, 84])]

# exact parametric-bootstrap null distribution of the clock likelihood ratio
def quick_LR(tt):
    nll, _ = make_nll(tt, main["vmeas"])
    full = min((optimize.minimize(lambda p: nll(p[0], p[1]), p0, method="L-BFGS-B",
                                  bounds=[(0, 60), (0, 60)]).fun for p0 in [(.3, .3), (1e-9, 1.)]))
    null = optimize.minimize_scalar(lambda u: nll(u, 0.0), bounds=(0, 60), method="bounded").fun
    return 2.0 * (null - full)
LR_null = np.array([quick_LR(T) for T in sim(SC0, 0.0, B_LR)])
p_LR_exact = float(((LR_null >= main["LR_clock"]).sum() + 1) / (B_LR + 1))

# ------------------------------------------------------ Dicke's ratio (1978)
def dicke(T):
    """Di = sum_{i=2}^n delta(i)^2 / sum_{i=2}^n (delta(i)-delta(i-1))^2, with
    delta the residual from the least-squares line through the n epochs."""
    T = T[None, :] if T.ndim == 1 else T
    n = T.shape[1]; A = np.column_stack([np.ones(n), np.arange(n, dtype=float)])
    d = T - (A @ np.linalg.lstsq(A, T.T, rcond=None)[0]).T
    return (d[:, 1:] ** 2).sum(1) / (np.diff(d, axis=1) ** 2).sum(1)

NM = len(ep_min)
di_obs = float(dicke(ep_min)[0])
di_ck, di_rw = dicke(sim(SC0, 0.0, B_SIM)), dicke(sim(main["sc"], main["sw"], B_SIM))
p_di_ck = float(((di_ck >= di_obs).sum() + 1) / (B_SIM + 1))
p_di_rw = float((2 * min((di_rw <= di_obs).sum(), (di_rw >= di_obs).sum()) + 1) / (B_SIM + 1))

# ------------------------------------------------- phase structure function
KMAX = 12
ks = np.arange(1, KMAX + 1, dtype=float)
def sfun(T, P):
    """Phase structure function, with the drift removed by the endpoint estimator
    -- which is the generalised-least-squares estimator of P under a random walk
    -- applied identically to observed and simulated records."""
    T = T[None, :] if T.ndim == 1 else T
    return np.array([np.mean((T[:, k:] - T[:, :-k] - k * P) ** 2, axis=1)
                     for k in range(1, KMAX + 1)]).T
wts = np.arange(len(ep_min) - 1, len(ep_min) - 1 - KMAX, -1).astype(float)
A = np.column_stack([np.ones(KMAX), ks]) * np.sqrt(wts)[:, None]
wls = lambda S: np.linalg.lstsq(A, S * np.sqrt(wts), rcond=None)[0]
S_obs = sfun(ep_min, (ep_min[-1] - ep_min[0]) / NC)[0]
S_rw = sfun(*(lambda T: (T, ((T[:, -1] - T[:, 0]) / NC)[:, None]))(sim(main["sc"], main["sw"], 4000)))
S_ck = sfun(*(lambda T: (T, ((T[:, -1] - T[:, 0]) / NC)[:, None]))(sim(SC0, 0.0, 4000)))
c_obs = wls(S_obs)
sl_rw = np.array([wls(s)[1] for s in S_rw]); sl_ck = np.array([wls(s)[1] for s in S_ck])
p_sf_ck = float(((sl_ck >= c_obs[1]).sum() + 1) / (len(sl_ck) + 1))

# ------------------------------------------ rigid ephemerides and prediction
n_idx = np.arange(len(ep_min), dtype=float)
res_free = ep_min - np.polyval(np.polyfit(n_idx, ep_min, 1), n_idx)
res_tid = ep_min - (P_TIDAL * n_idx + np.mean(ep_min - P_TIDAL * n_idx))
pred = float(ep_min[-1] + main["P"])
pred_se = float(np.sqrt(main["sw"] ** 2 + main["sc"] ** 2 + main["P_se"] ** 2
                        + sig_min[-1] ** 2 + QV))
ncoh = float((main["P"] / main["sw"]) ** 2)
tot_sig = np.sqrt(sig_min ** 2 + QV)
L_obs = np.diff(ep_min)

res = {
    "n_months_total": int(len(t)), "n_months_missing_ssn": int(np.isnan(x).sum()),
   "n_months_with_stated_error": int(has_err.sum()),
   "n_months_provisional": int((~df["definitive"].to_numpy(bool)).sum()),
   "record_start": float(t[0]), "record_end": float(t[-1]),
   "median_monthly_se": float(np.median(se[has_err])),
   "median_monthly_se_1818_1849": float(np.median(se[has_err & (yr <= 1849)])),
   "median_monthly_se_2000_2026": float(np.median(se[has_err & (yr >= 2000)])),
   "n_minima": int(len(ep_min)), "n_maxima": int(len(ep_max)),
   "first_minimum": float(ep_min[0]), "last_minimum": float(ep_min[-1]),
   "epoch_sigma_median": float(np.median(tot_sig)), "epoch_sigma_max": float(tot_sig.max()),
   "epoch_sigma_max_year": float(ep_min[int(np.argmax(tot_sig))]),
   "n_epochs_unmoved_by_noise": int((sig_min < 1e-6).sum()),
   "epoch_quantisation_sigma": float(np.sqrt(QV)), "n_cycles": main["ncyc"], "P": main["P"],
   "P_se": main["P_se"], "mean_cycle_length": main["meanL"], "sd_cycle_length": main["sdL"],
   "shortest_cycle": float(L_obs.min()), "longest_cycle": float(L_obs.max()),
   "sigma_w": main["sw"], "sigma_w_lo68": main["sw_lo"], "sigma_w_hi68": main["sw_hi"],
   "sigma_c": main["sc"], "sigma_c_ul95": main["sc_ul95"], "sigma_c_under_clock": SC0,
   "LR_clock": main["LR_clock"], "p_clock_asymptotic": main["p_clock"],
   "p_clock_bootstrap": p_LR_exact, "n_bootstrap_LR": B_LR,
   "significance_sigma_clock": float(stats.norm.isf(max(main["p_clock"], 1e-300))),
   "LR_randomwalk": main["LR_rw"], "p_randomwalk": main["p_rw"], "r1": main["r1"],
   "r1_lo68": r1_lo, "r1_hi68": r1_hi, "r1_clock_expectation": -0.5,
   "r1_clock_sim_mean": float(r1_ck.mean()), "r1_clock_sim_sd": float(r1_ck.std(ddof=1)),
   "p_r1_vs_clock": p_r1_ck, "r1_rw_sim_mean": float(r1_rw.mean()),
   "r1_rw_sim_sd": float(r1_rw.std(ddof=1)), "p_r1_vs_randomwalk": p_r1_rw,
   "dicke_ratio": di_obs, "dicke_n": int(NM),
   "dicke_expect_clock": float((NM ** 2 - 1) / (2 * (NM ** 2 + 2 * NM + 3))),
   "dicke_expect_randomwalk": float((NM + 2) * (NM ** 2 - 1) / (3 * (5 * NM ** 2 + 6 * NM - 3))),
   "dicke_clock_sim_mean": float(di_ck.mean()), "dicke_clock_sim_sd": float(di_ck.std(ddof=1)),
   "dicke_rw_sim_mean": float(di_rw.mean()), "dicke_rw_sim_sd": float(di_rw.std(ddof=1)),
   "p_dicke_vs_clock": p_di_ck, "p_dicke_vs_randomwalk": p_di_rw,
   "dicke_ratio_maxima": float(dicke(ep_max)[0]), "sf_kmax": KMAX, "sf_slope": float(c_obs[1]),
   "sf_slope_sim_sd": float(sl_rw.std(ddof=1)), "sf_intercept": float(c_obs[0]),
   "sf_clock_sim_slope_mean": float(sl_ck.mean()),
   "sf_clock_sim_slope_sd": float(sl_ck.std(ddof=1)),
   "sf_rw_sim_slope_mean": float(sl_rw.mean()), "p_sf_vs_clock": p_sf_ck,
   "n_bootstrap_LR_exceeding": int((LR_null >= main["LR_clock"]).sum()),
   "endpoint_period": float((ep_min[-1] - ep_min[0]) / NC), "sigma_w_no_meas_err": no_err["sw"],
   "meas_var_fraction_of_sigma_w2": float(np.mean(main["vmeas"]) / main["sw"] ** 2),
   "maxima_n_cycles": maxima["ncyc"], "maxima_P": maxima["P"], "maxima_P_se": maxima["P_se"],
   "maxima_sigma_w": maxima["sw"], "maxima_sigma_c": maxima["sc"],
   "maxima_LR_clock": maxima["LR_clock"], "maxima_p_clock": maxima["p_clock"],
   "post1818_n_cycles": post["ncyc"], "post1818_P": post["P"], "post1818_P_se": post["P_se"],
   "post1818_sigma_w": post["sw"], "post1818_sigma_c": post["sc"],
   "post1818_LR_clock": post["LR_clock"], "post1818_p_clock": post["p_clock"],
   "smooth24_P": sm24["P"], "smooth24_sigma_w": sm24["sw"], "smooth24_sigma_c": sm24["sc"],
   "smooth24_LR_clock": sm24["LR_clock"], "smooth24_p_clock": sm24["p_clock"],
   "LR_period_drift": float(LR_drift),
   "p_period_drift": float(stats.chi2.sf(max(LR_drift, 0.0), 1)),
   "clock_resid_rms": float(np.sqrt(np.mean(res_free ** 2))),
   "clock_resid_max_abs": float(np.max(np.abs(res_free))),
   "clock_resid_max_year": float(ep_min[int(np.argmax(np.abs(res_free)))]),
   "tidal_period": P_TIDAL, "tidal_resid_rms": float(np.sqrt(np.mean(res_tid ** 2))),
   "tidal_resid_max_abs": float(np.max(np.abs(res_tid))),
   "P_minus_tidal_sigma": float((main["P"] - P_TIDAL) / main["P_se"]),
   "pred_cycle26_minimum": pred, "pred_cycle26_minimum_se": pred_se,
   "phase_drift_per_cycle_deg": float(360.0 * main["sw"] / main["P"]), "coherence_cycles": ncoh,
   "coherence_years": float(ncoh * main["P"]),
   "phase_diffusion_cycles2_per_year": float(main["sw"] ** 2 / main["P"] ** 3),
   "hypothesis_supported": bool(p_LR_exact < 0.05 and main["sw_lo"] > 0 and post["p_clock"] < 0.05 and maxima["p_clock"] < 0.05)
}
for k in (1, 2, 4, 8, 12):
    res[f"sf_obs_k{k}"] = float(S_obs[k - 1])
    res[f"sf_rw_k{k}"] = float(S_rw[:, k - 1].mean())
    res[f"sf_clock_k{k}"] = float(S_ck[:, k - 1].mean())
for k, v in res.items():
    if isinstance(v, float) and not np.isfinite(v):
        raise SystemExit("non-finite: " + k)
print("minima :", np.round(ep_min, 3).tolist())
print("sigma  :", np.round(tot_sig, 4).tolist())
print("lengths:", np.round(L_obs, 3).tolist())
print("S(k)   :", np.round(S_obs, 3).tolist())
print("RESULT_JSON: " + json.dumps(res))

Results and analysis

Key numbers

Caveats

• Only 24 cycle intervals exist in the entire directly measured record, so the 68 per cent interval on sigma_w spans 1.01-1.34 yr and the upper limit on any clock component is correspondingly weak; the study bounds a chronometer rather than excluding one of arbitrarily small amplitude.
• The test distinguishes a memoryless random walk from a rigid clock. A slow deterministic modulation of the period, of Gleissberg type, would also make the phase structure function rise, and with 24 cycles it cannot be cleanly separated from genuine diffusion; the null of no linear period drift (p = 0.62) rules out only the simplest such alternative.
• Sunspot minima are a proxy for the phase of the underlying magnetic cycle, and near minimum the smoothed series is flat, so the marker is intrinsically noisy. The maxima cross-check gives the same verdict but a larger sigma_w (1.958 yr), showing that the numerical value of the diffusion constant is marker-dependent even though the qualitative conclusion is not.
• Only random measurement error is propagated. Systematic errors in the early SILSO version 2.0 record, where a single observer often provides all monthly values, could shift epochs coherently in a way the Monte Carlo does not capture; the 1823-onward subsample is the partial control on this and gives a smaller sigma_w (0.855 yr).
• The epoch definition quantises to the month and the flat 1810 minimum required a plateau-averaging rule, so individual epochs carry a definitional as well as a statistical uncertainty; the 24-month-smoothing variant is the check on this and changes sigma_w by about 5 per cent.
• The parametric-bootstrap p-value is limited by the 2000 replicates to p < 5e-4; the quoted 5.2 sigma comes from the asymptotic boundary-corrected mixture distribution, which the bootstrap corroborates but cannot resolve.
• The last five months of the record are provisional and the September 2026 file will be revised, but the most recent epoch used is the 2019.96 minimum, so no result depends on provisional values.

≈ 2552 words · AI-drafted manuscript · not peer-reviewed

A compact, well-executed Letter that settles a genuinely contested question with the right statistical tool and honest error propagation. The phase model of Eqs. (1)-(2) nests the rigid clock and the random walk, and on the 25 sunspot-cycle minima of the SILSO v2.0 record (3332 monthly means, 1749.042-2026.623, WDC-SILSO, retrieved 2026-09-08) the clock limit is rejected with 2 dlnL = 27.04 while sigma_c sits on its boundary. Every manuscript number I checked - sigma_w = 1.155 (1.006-1.344), sigma_c < 0.156, rho_1 = +0.103, Di(25) = 2.171 against 0.460 and 1.716, S(1) = 1.333 and S(12) = 10.191 yr^2, the maxima/1823-onward/24-month variants at 1.958/0.855/1.092 yr, the drift null at p = 0.62, and the 2031.0 +/- 1.2 forecast - matches the executed script and its recorded output, the sample (25 minima, 24 intervals, 24 maxima; 19 epochs from 1823) matches the code, and the data are correctly and never misleadingly described as the real SILSO series. The Letter is within length and follows the venue's conventions. What holds it back from acceptance as it stands is the treatment of uncertainty on the central number: the analysis-choice spread in sigma_w (0.855-1.958 yr) is larger than the quoted statistical interval and is not carried through to the diffusion constant or the coherence horizon, the 5.2 sigma is asymptotic and the bootstrap resolves only to 3.5 sigma, and the bound on sigma_c is stated more absolutely than the model supports. These are matters of framing and disclosure, not of re-analysis, and the qualitative conclusion - no chronometer, a diffusing phase - survives every variant the authors ran.