Changelog

This page renders the project’s root CHANGELOG.md directly, so the changelog has a single source.

All notable changes to Magνs are documented in this file.

The format follows Keep a Changelog, and the project uses Semantic Versioning.

[Unreleased]

Added

• The paper, arXiv:2610.07159. The README and the docs landing page carry an arXiv badge. “How to cite”, in the README and on the docs’ citing page, asks for the paper and for the version of the software used, and the citing page gives the paper’s BibTeX entry as INSPIRE holds it. CITATION.cff names the paper as its preferred-citation, so GitHub’s “Cite this repository” offers it.

• A Zenodo DOI, 10.5281/zenodo.23160508. This concept DOI resolves to the latest version; version 1.2.0 has its own, 10.5281/zenodo.23160509. The README and the docs landing page carry a DOI badge, CITATION.cff lists both DOIs, and the software’s BibTeX entry on the docs’ citing page gives the concept DOI.

[1.2.0] - 2026-10-05

Fixed

• Notebook 28’s long-range solar figure read the probability out inside the Sun. The L_e - L_mu sweep ended at the edge of the solar table, 0.98 R_sun, where the potential of a mediator of range R_sun is still about 1e3 V_CC. Read out in the eigenbasis there, the curve turned up above about 9 MeV and crossed the standard one, an effect no detector on Earth sees: the potential fades within a few solar radii and the exit is adiabatic. The sweep now ends at 20 R_sun, where the readout is the vacuum one. It also starts at 0.05 R_sun, where 8B neutrinos are made, instead of the center, and uses one coupling for both mediator ranges, so that the two curves differ only through the range. The figure’s cache entry is recomputed.

• Notebook 28 draws the 8B flux at Earth under the long-range force, on B16-GS98. The long-range figure moves from BS2005-AGS,OP, whose table stops at 0.98 R_sun, to B16-GS98, which reaches the surface. A new figure, solar_8b_flux.pdf, gives the nu_e flux at Earth from 8B decay for the standard case and both mediator ranges: the probability averaged over where 8B is made, times the spectrum of Winter et al. (2006) and the B16-GS98 total flux. A dotted curve triples g’^2 at 1/m = R_sun. The production profile is rebuilt from the model’s structure, since the published one is no longer online; docs/dev/measurements/solar_8b_flux/ holds the inputs, the script, and its check against Bahcall’s published BS05 distribution.

• The plotting functions’ axis and tick labels were 10 points outside the notebooks. magnus.plotting left text sizes to the matplotlibrc of wherever it ran, so in the documentation and in users’ scripts the 18-by-9-inch figures had tiny labels. The module now brings its own sizes, magnus.plotting.HOUSE_RC (25-point axis labels, 23-point tick labels, ticks pointing in on all four sides), and reads no matplotlibrc; a size the caller has changed is kept. Two layouts were adjusted to fit the larger text: plot_biprobability now calls tight_layout by default (its ordinate label was cut off), and a probability panel that plot_probability_with_profile labels itself shows the probability symbol only (the long label ran into the density label).

• Warnings that say “Shown once per session” are shown once per session (issue #205). Python’s default filter remembers a message per call site and forgets it whenever any code changes the filters, which the package itself does during a calculation; one run of the paper’s cavity example printed the same MagnusConvergenceWarning four times. These messages (MagnusConvergenceWarning, ScalarHamiltonianWarning, ToleranceNotAchievedWarning and its subclasses) now go through one helper that remembers what it has shown. A filter the caller sets still decides: 'always' shows every one, 'error' raises, 'ignore' hides.

• The windowed average counts the propagations that warned, not the warnings (issue #193). One propagation can raise the same class twice, so the summary read “Raised by 10 of the 5 propagations”. It now counts each propagation once.

• Documentation audit (issue #193). Printed values brought up to date with #168, broken and incomplete listings fixed, contradictions between pages and with the code removed, and development history moved out of the user pages.

• average=True no longer returns zero at small phase on a constant Hamiltonian (issue #163). In vacuum or at constant density, the phase average was computed and then discarded for the zero-phase limit wherever the two differed by less than 1e-4 in absolute terms, whatever rtol and atol: at 1 GeV over 1 to 3 km, P(numu -> nue) = 4.7e-7 came back as 1e-32, with no warning. The comparison is now relative as well as absolute, and follows the caller’s tolerance on every route, as it already did on a profile. Results that were right are unchanged, bit for bit; the averaged values in the docs and the paper do not move.

• A mass density becomes an electron density as ρ N_A Y_e (issue #168). The package divided ρ Y_e by the mean free-nucleon mass, (m_p + r m_n)/(1 + r), which ignores nuclear binding: matter weighs the atomic mass unit per nucleon, 931.494 MeV, not 938.9 MeV, so every electron density converted from g/cm³ was 0.8% low. It also converted grams with a rounded e/c², 1.783e-33 g/eV, another 1.9e-4. The Earth’s Y_e values were already electrons per atomic mass unit (iron’s 26/55.845, seawater’s 10/18.015), so for them ρ N_A Y_e is exact; nuSQuIDS, GLoBES, Prob3++, NuFast, OscProb and nuCraft all use it. Now matter.num_density_e_func divides by the new globaldefs.ATOMIC_MASS_UNIT, CONV_EV_TO_G is exact, and the neutron-to-proton ratio, which entered the conversion only through the nucleon mass, no longer does (it still sets the sterile states’ matter entry). The tabulated solar models count electrons, and the neutron-to-proton ratio, from the atomic masses of ¹H and ⁴He, the rest counted as helium: within 0.07% of BS05(OP)’s full composition, where ρ(1 + X)/(2 m_N) was 0.1 to 0.4% low and the textbook ρ N_A (1 + X)/2 is 0.4 to 0.8% high. Every result computed from a mass density changes: through the Earth’s core by up to 6.5e-2 in probability at 1-20 GeV, through the mantle 1e-2, at DUNE’s constant density 5e-4; averaged solar probabilities by up to 4.5e-4, and the exponential solar fit, defined as an electron density with N_A, not at all. The supernova-shock model of the notebooks and the benchmark PREM chord keep their normalization, as model definitions their frozen references were built on.

• The LIV term changes sign for antineutrinos when n_liv is even (issue #162). For antineutrinos the Lorentz-violating operator is transposed in flavor space and its CPT-odd part changes sign; an operator of dimension n_liv + 3 is CPT-odd when that is odd (Kostelecky & Mewes, Phys. Rev. D 85, 096005 (2012), Eq. 78). nubar=True only conjugated the term, whatever n_liv: right for odd n_liv, the wrong sign for even n_liv, the default n_liv=0 included. At two flavors nubar did not reach the LIV builder at all. Now every LIV builder and wrapper applies the sign, hamiltonian_2nu_liv_energy_independent takes nubar, and osc_prob_2nu_vacuum_liv declares it, so the CLI’s --nubar, which it used to drop there without a word, now reaches it. Antineutrino results at even n_liv change: on the paper’s Fig. 13 chord ($\cos\theta_z = -0.9$, $b_3 = 5.4\cdot10^{-14}$ eV, n_liv=0) the $\bar\nu_\mu$ survival probability at 1.46 GeV goes from 0.536 to 0.353. Neutrinos and odd n_liv are unchanged, and no notebook or paper figure computes LIV for antineutrinos. The docs no longer call every LIV term CPT-odd.

• An energy scan whose slab count starts at the cap says it checked nothing (issue #184). Across the Sun at rtol=atol of 1e-4 or tighter, every energy’s slab count starts at max_n_slabs, so the energy-batched ladder computes one level with nothing to compare it with; the warning said only that a refinement cap was reached. It now says the answers are unverified, that across the Sun about 1e-4 is the tightest tolerance a scan can verify, and to compare with strategy='hybrid', loosen the tolerance, or raise max_n_slabs. Over 0.5-20 MeV the returned level was up to 2.0e-3 off on the default profile, where the hybrid was within 6.3e-5. No answer changes. diagnostics.rst documents it, and documents that the hybrid’s certified=True is an agreement between levels, like the ladder’s: on B16-GS98 at 0.7 MeV it certified 1e-6 on an answer 1.8e-5 off (issue #187).

• cross_check_strategies says which engines did not reach the tolerance (issue #166). On the Sun at 1 MeV it reported a 1.9e-3 spread between the certified hybrid and the Magnus ladder, and named the pair. Against a 2e6-slab reference the hybrid is 3.8e-5 off and the ladder, stopped at its 20 000-slab cap on a single level, 1.9e-3; the ladder had warned ToleranceNotAchievedWarning, the hybrid nothing. New keys: 'unverified', the engines that raised ToleranceNotAchievedWarning or HybridCertificationWarning, and 'max_spread_verified' with its pair, the spread across families among the others. The existing keys and values are unchanged, and no engine runs that did not before.

• A tight tolerance on the Sun no longer returns a worse answer than the default (issue #167). At 1 MeV across the Sun, rtol=atol=1e-8 is beyond what the adiabatic hybrid certifies, so strategy='auto' declined it and the Magnus ladder answered instead; that ladder needs about 2e6 slabs, stopped at max_n_slabs = 20 000 on a single level, and returned it 4.0e-3 off, where the default tolerance gives 8.8e-6. (The issue read this as 'auto' choosing the ladder up front; it tries the hybrid first.) The hybrid’s uncertified answer is now kept when it declines, and a point the per-point ladder then fails to converge gets it back when the ladder has less to show for its own: no two levels to compare, or two further apart than the hybrid’s error estimate (GAMMA_TO_ERROR times the largest non-adiabaticity along the path). ToleranceNotAchievedWarning says so, and strategy_info names the hybrid, uncertified, with the estimate. Against a 2e6-slab reference at 1e-8: 4.0e-3 → 1.7e-6 at 1 MeV, 7.9e-4 → 1.9e-5 at 10 MeV, 1.9e-3 → 7.4e-8 on B16-GS98, 9.3e-4 → 1.1e-6 with NSI, 5.5e-3 → 1.4e-6 at four flavors. Calls that converge are untouched: 139 of 140 benchmark calls are bitwise identical, and the 140th, the 5-flavor Sun at 8 MeV and the default tolerance, moved from 2.9e-4 to 3.5e-6 off a converged reference. Still open: a scan of several energies at such a tolerance is answered by the energy-batched engine, which this does not reach (4.0e-3 at 1 and 10 MeV).

• Undeclared jumps and narrow spikes no longer end the refinement ladder on a wrong answer (issue #161). One point at a time, a slab straddling a density jump nobody declared converges as 1/n, and the ladder stopped wherever two levels happened to agree: the castle wall of the issue was off by 4.3e-2, and 24 of 40 baselines by more than the requested 1e-3, with only MagnusConvergenceWarning. osc_prob_energy_baseline now watches the ladder’s own samples: where the largest change between nodes does not shrink from one level to the next, it locates the jumps on the hybrid strategy’s 200-point probe grid, declares them, recomputes, and warns UnmarkedDiscontinuityWarning naming them (4.1e-2 → 4.5e-13). A narrow smooth spike was stepped over by every level instead: the 5 km spike of #154 was off by 5.4e-3 with no warning at all, and 17 of 36 spikes of 2-50 km by more than 1e-3. strategy='auto', which probes the profile anyway, now starts the ladder at the probe spacing where the probe flags a sharp feature (AUTO_SHARP_FEATURE_SLABS_PER_PROBE): 17 → 3, the three a 2 km spike between two probes, which no grid of that spacing can see. Of 140 public-API calls benchmarked against 1.1, 126 are bit-identical and the other 14 moved toward the reference. Timing: no change on any path with declared breakpoints (every Earth wrapper), on the Sun under 'auto', or in vacuum and constant density; +2 % to +5 % (about 0.05 ms, 0.7 ms on a 20 ms solar ladder) on a single point that reaches the ladder without breakpoints, the cost of watching its samples. The issue’s first proposal, treating an agreement to rounding as unproven, was measured and not adopted: it caught 19 of the 41 wrong answers and fired on 6 of 27 right ones.

• average=True on the 5-flavor Sun at its default sterile parameters returns (issue #148). There the two sterile states are decoupled and degenerate along the whole path; the adiabatic averaging route read the pair as a crossing spanning the Sun, laid 338 193 energy nodes across it, and never returned (killed after 90 s; a 40-energy scan ran past 5 minutes). A flavor that nothing couples to is an eigenstate at every position and never oscillates, so when one is degenerate with another level along the whole path it is now averaged out: probability 1 of staying itself, 0 to or from anything else, and the rest averaged as before. The 5-flavor Sun at its defaults, with or without NSI, on the default profile or BP04, now returns in about 0.01 s the 3-flavor value to 2e-15. Calls without such a pair are unchanged, bit for bit – the 4-flavor Sun at its defaults included, whose sterile state crosses an active level at a point, which the windows handle (#59). As a backstop, the grid route of the phase average now stops at avgprob.PHASE_AVERAGE_MAX_GRID_NODES = 20 000 nodes, 250 times the most any test takes, and the call returns the decohered limit with PhaseAveragingWarning.

• Calls from several threads no longer interfere (issue #153). Four pieces of per-call state were module globals: the two slab-norm collectors of magnus.magnus and the engine trace and disabled set of magnus.oscprob. Concurrent batched Earth energy scans overwrote each other’s collector – over 100 threaded calls, 14 to 20 raised ValueError: need at least one array to concatenate and 7 more returned different numbers – and strategy='magnus' in one thread disabled the other engines for every thread while it ran (30 of 40 mixed-strategy calls differed from serial). They are contextvars.ContextVars now, so each thread has its own; concurrent calls match the serial answer bit for bit. performance.rst says so, and that n_jobs uses processes.

• The sdist’s own tests pass against the installed package (issue #164 §1). Thirty-eight tests check the repository rather than the library – the documentation, the notebooks, the paper’s assets, the CI workflows – and the sdist ships none of those files, so run from an unpacked sdist they failed (20) or errored (18). They are marked checkout_only and skip there with their reason. tests/conftest.py, which holds that rule, had never been in the sdist either; MANIFEST.in now includes it. A new CI job, sdist tests, builds the sdist and the wheel, installs the wheel into a clean environment and runs the sdist’s tests with src/ removed. In a checkout nothing is skipped.

• average_spread and average_initial_state are refused without average=True (issue #160). They shape the phase average and did nothing without it; they were documented as ignored and are now refused by name, like average_n_samples. osc_prob_energy_baseline went on ignoring average_n_samples without average too, since its direct route reaches the averaging code only with average=True; it refuses all three now. Left at their default, None, they pass as before.

• Nothing calls numba optional any more (issue #164 §2). numba and matplotlib have been base dependencies since the fast and plot extras were emptied, but HAVE_NUMBA’s docstring still called numba “an optional dependency (pip install magnuspy[fast])”, and the error for expm_backend='numba' without numba pointed at that extra. Both now say numba is a dependency and name pip install numba. The empty plot and fast extras are kept, so pip install 'magnuspy[plot]' from older instructions still resolves.

• The plotting tests pass at the declared matplotlib floor, 3.5 (issue #164 §3). Three tests read API from later releases (Figure.get_supylabel, the figure.labelsize rcParam, the legend’s _ncols), and one relied on 3.6’s legend keeping more columns than entries, so they failed at matplotlib 3.5.0 while the library itself worked. They now check the same behavior through what 3.5 has; the runtime floor is unchanged.

• Averaging an energy-independent Hamiltonian says why no average applies (issue #144 §3). A matrix, or a function of position alone, has no energy dependence for a spread to act on, so a pair of levels neither decohered nor coherent stays that way; the warning said only that such a pair existed. It now says the Hamiltonian does not depend on energy, that no spread can decohere the pair, and that average=False is the meaningful quantity.

• The energy-window average raises one warning per class, not one per sample (issue #144 §2). On a profile with declared discontinuities, average=True propagates once per energy sample, and each propagation warned: 41 MagnusConvergenceWarning and 41 ToleranceNotAchievedWarning records for one call. They are collected and raised once per class, saying how many of the propagations raised it; filters set by the caller still apply.

• “Shown once per session” is said only where it is true (issue #144 §1). Python’s default filter shows each distinct message once per call site, so a warning whose text reports a count, a window, an error or a slab number shows again for each new value. The phrase is gone from the nine such warnings (the averaging, hybrid-certification, energy-batched and cumulative-scan ones, and MagnusHighOrderCostWarning) and kept on the fixed-text ones and on those whose text is one of a few fixed phrases by design. diagnostics.rst says the same.

• Compute-mode plots name the plotting routine in a wrapper’s error (issue #160 §12). A mistyped physics argument was refused naming only oscprob.osc_prob_3nu_earth, a function the caller never called; the message now starts with the plotting routine and keeps the wrapper’s own text.

• A computed oscillogram refuses log10_energy above 19, the Planck scale in GeV (issue #160 §12). [30, 31] was computed as if meaningful.

• plot_probability_vs_energy warns when the energies look like eV under a larger unit (issue #160 §12). Energies left in eV, the package’s unit, were drawn under the default GeV label a billion times too large. EnergyUnitWarning fires when the values reach 1e6 and would put the axis above 100 PeV.

• A curve with nothing above 0 is refused on a log y axis (issue #160 §12). plot_curves and plot_curves_stacked drew an empty axis for it, without a word. plot_curves_stacked now checks each panel’s curves as plot_curves does.

• A scalar-only rho_func is warned about by name, once (issues #144 §4, #160 §7). The engine’s ScalarHamiltonianWarning told the caller to rewrite H_func, a function they never wrote. The scenario functions now probe rho_func at three points; if it cannot take an array, the warning names rho_func and the NumPy fix, and it is evaluated position by position so the Hamiltonian stays array-capable and the engine does not warn again. Per-point results are bit-identical (14 of 14 cases across the three scenario functions and both strategies). An energy scan at one baseline now stays on the batched engine instead of falling back to one ladder per energy: it moves within the tolerance (4.4e-15 at the default, 6.5e-8 at rtol=1e-6 on a 40-energy scan) and runs about twice as fast at the default tolerance.

• adiabatic.hybrid_propagator refuses a Hamiltonian that turns NaN along the path (issue #160 §7). Only the first sample was checked, so a NaN past it came back as a NaN propagator, silently.

• A custom Sun or Earth Hamiltonian is checked at its first sample (issue #160 §7). osc_prob_sun ran an H_func returning a list; it is now refused by name, as it already was by osc_prob_energy_baseline, and a non-square, non-finite or non-Hermitian first sample is refused too. One evaluation, at L0 and the first energy.

• osc_prob_energy_baseline refuses an H_func taking no argument, or returning an object array, by name (issue #160 §7). The first failed on its first call as a TypeError naming a lambda, the second as NumPy’s “ufunc ‘isfinite’ not supported”.

• t_breakpoints of which none lies on the path are refused (issue #160 §6). Breakpoints all before L0 or past the longest baseline mark nothing – a sign or a unit gone wrong – and were accepted and ignored. Some off the path stay accepted: a whole profile’s breakpoints reused for a shorter path, as solar_models.rst does with a table’s rows, still mark the jumps on it. The Earth wrappers now drop the PREM crossings past the end of a partial path themselves; results are bit-identical (30 of 30 Earth cases, full and partial paths, scans and strategies).

• average=True with cumulative=True is refused (issue #160 §5). The phase average takes no baseline scan, so the explicit request for one was accepted and ignored.

• The n_jobs docs say the worker pool outlives the call (issue #160 §5). joblib’s loky workers stay alive until idle for 300 s, which looked like leaked processes.

• n_jobs=None is refused by name (issue #160 §5). It was let through the check and then compared with an integer, as “’<’ not supported between instances of ‘int’ and ‘NoneType’”; on single-point calls it was accepted. n_jobs is -1 or a positive integer.

• osc_prob documents what happens to a count above its cap (issue #160 §5). With a tolerance, an n_slabs floor above max_n_slabs is clipped to it and the points per slab are clipped at max_n_tpts_per_slab, with ToleranceNotAchievedWarning if a cap stops the refinement; on a fixed grid (rtol=atol=None) both counts are used as given. This was the behavior already, undocumented; it is unchanged.

• A filename_log that cannot be written is refused by name, before any work (issue #160 §4). With save_log=True, a missing directory or a path naming a directory surfaced as a raw FileNotFoundError or IsADirectoryError from deep inside osc_prob.

• default_osc_params_set_name=None is refused by name (issue #160 §4). It failed as “can only concatenate str (not NoneType) to str”; the name of a parameter set is a string.

• osc_prob checks new_recursion_limit (issue #160 §4). The scenario functions already refused a non-positive or non-integer value; osc_prob accepted 0, -1, 2.5 and 'a'. It is still otherwise ignored.

• A points-per-slab setting passed with integration_method='gl' warns that it does nothing (issue #160 §5). n_tpts_per_slab, min_n_tpts_per_slab, max_n_tpts_per_slab and growth_factor_n_tpts_per_slab are used by 'trapezoid' and 'simpson' only; 'gl' evaluates the Hamiltonian at its own nodes and overrode them silently. A value the caller passes now raises the new IgnoredQuadratureSettingWarning, once per call, at the caller’s line; defaults forwarded between layers never do. The result is unchanged. The docstrings, methodology.rst and diagnostics.rst say so, and the notebooks no longer pass these settings with 'gl' (68 calls; outputs unchanged).

• An odd magnus_exp_order with integration_method='gl' is refused (issue #160 §5). The Gauss-Legendre schemes have even orders only, and an odd order ran the next even one bit for bit (checked on main at orders 1, 3, 5 and 7 through every engine): magnus_exp_order=3 was order 4 under another name. The error names the even order. Odd orders still run on 'trapezoid' and 'simpson'. A constant Hamiltonian’s single-slab shortcut now runs at order 2 instead of 1, with results bit-identical to before. The notebooks that passed odd orders now pass the even order they were running, so their outputs are unchanged; notebook 24’s order table shows 2, 4 and 6.

• Flags given as None, and validate_input given as anything but a bool, are refused (issue #160 §4). The entry checks skipped every None, so nubar=None, average=None, cumulative=None and seven more flags read as False; validate_input was truth-tested before anything was checked, so None turned validation off and an array raised NumPy’s “truth value is ambiguous”.

• L0 on a vacuum or Earth entry point is refused by name (issue #160 §1). None of them declares it, and one passed anyway collided with the start the entry point sets, as a TypeError naming an internal function; the vacuum LIV wrappers used it as the start.

• The vacuum wrappers refuse t_breakpoints (issue #160 §1). The vacuum Hamiltonian has no discontinuity to mark, so breakpoints were accepted and ignored. The refinement keywords (rtol, atol, the slab and order controls) stay accepted, since the command line and shared calls pass one set to every wrapper, and the vacuum docstrings now say they have no effect there.

• average_n_samples without average=True is refused (issue #160 §1). It was accepted and ignored; it was already refused on the averaging routes that do not sample. average_spread and average_initial_state stay documented as ignored without average.

• The scenario functions refuse an h_vac_energy_indep they would ignore (issue #160 §1). Up to five flavors the vacuum Hamiltonian is built from osc_params, so a matrix passed there was accepted and ignored; osc_prob_energy_baseline is the way to propagate one’s own.

• The Earth wrappers refuse a neutron-to-proton ratio they ignore, and check a callable one (issue #160 §1). The density derives its ratio from :math:Y_e, so the caller’s enters only the sterile projector: at two and three flavors a ratio, scalar or callable, was accepted and ignored, and is refused now. Above that, a callable ratio returning NaN gave NaN probabilities; its value at the start of the path is checked, as rho_func’s is.

• The wrappers broadcast a single-entry energy or L against a longer one (issue #160 §1). osc_prob_energy_baseline always did; the osc_prob_* wrappers refused [1e9] against three baselines as a length mismatch. Unequal lengths above one are still refused.

• The Sun wrappers refuse a negative L0 and arguments with no effect (issue #160 §1). L0 is the radius where the path starts; a negative one was accepted. electron_fraction was accepted and ignored by every Sun wrapper, which uses the solar model’s electron density directly, and ratio_number_neutrons_to_protons was accepted and ignored at two and three flavors. All three are refused by name.

• A NaN l_scale or rho_central on an exponential-density wrapper is named as such (issue #160 §1). It was refused only once the profile returned NaN, as “rho_func must be finite”, naming an argument the caller never passed.

• earth.earth_radial_distance_from_depth checks the position l (issue #160 §3). A NaN, negative or True position was accepted (NaN returned NaN). It is refused now, at the cost of one array comparison, which also catches -inf.

• earth.dms_to_decimal refuses out-of-range degrees and ambiguous signs (issue #160 §3). Degrees of 500 were converted, and (10, -1, 0) was read as 10°1′ south. The degrees now lie in [-360, 360], and a negative part after a positive one is refused, as the chord helpers already did: the sign goes on the first nonzero part.

• The per-node matter helpers check their fractions (issue #160 §11). earth.electron_fraction_func_prem, earth.neutron_to_proton_ratio_from_electron_fraction and matter.num_density_e_func accepted negative, NaN, infinite and True electron fractions (or a negative neutron/proton ratio, or a density that is not a function), and computed with them. A plain float in range costs one comparison; anything else is checked in full when the call comes from outside the package, where the wrappers have not already checked it. Earth calls time as before.

• MagnusConvergenceWarning measures the traceless part of Ω (issue #155 §1). The trace is a global phase: it commutes with every term of the Magnus series and has no bearing on its convergence, but it was counted, so a custom Hamiltonian’s first coarse slab was flagged (‖Ω‖ = 3.19, 1.94 without the trace) at every tolerance, and rtol=1e-8 did not silence it. Over the 160 valid configurations of warn_fp.py the warning now fires 19 times instead of 39, with 7 true positives instead of 8 and 12 false positives instead of 31; every answer outside the tolerance still carries a warning. The message now names n_slabs and min_n_slabs as what silences it; a smaller rtol refines, but the coarse first level is still reported. The refinement gate of the batched scan is unchanged.

• MixingAngleConventionWarning judges by the active angles, and says how to silence it (issue #160 §11). With angles='deg', sines in the active slots went unnoticed when sterile angles of a few degrees were also given, since the check looked at the largest of all angles. It now looks at θ₁₂, θ₁₃, θ₂₃ (or θ at two flavors), and its text says that genuinely tiny angles set it off and how to filter it.

• An angle beyond ±90° under angles='rad' or 'deg' is refused (issue #160 §11). The cosine is taken as +sqrt(1 - sin^2), so θ₁₂ = 2.0 rad silently gave the rotation of π − 2.0, 0.0175 off in probability. The message names the angle and the in-quadrant angle with the same sine.

• A computed profile plot labels its probability panel (issue #146 §2). plot_probability_with_profile in compute mode left the probability panel’s ordinate empty; it now gets the channel’s label, as plot_probability_vs_energy gives, unless panel_ylabels or shared_ylabel is passed.

• An oscillogram computed for antineutrinos is labeled as one (issue #145 §1). plot_oscillogram(..., wrapper_kw=dict(nubar=True)) computed the antineutrino grid but labeled its color bar and annotation with the neutrino channel.

• Small pseudo-Dirac splittings keep their precision (issue #165 §1). hamiltonian_pseudo_dirac_vacuum_energy_independent formed each splitting-dependent entry as (m2_j + delta) - m2_j at the scale of m2_j, so delta = 1e-18 eV^2 kept one or two digits (0.24 off at 1e-19 on the third state). It now builds the Hamiltonian as a base term, written down directly with nothing to cancel, plus the splitting term, carried at full precision: delta is recovered to 2e-16 on every state from 1e-19 to 1e-15. At ordinary splittings the result moves by rounding only (4e-19 on entries of 1e-3), and it is Hermitian to the last bit.

• The docs say that a parallel scan agrees with the serial one to within the tolerance, not bit for bit (issue #166 §1). n_jobs > 1 sends a scan to the per-point path, one refinement ladder per point, where the batched engine runs one for all of them: up to 3.7e-5 apart on a 40-energy Earth chord at the default 1e-3. Stated in the n_jobs docstring of osc_prob_energy_baseline and on the performance page.

• n_jobs > 1 no longer costs a second or more on a small scan (issue #155 §3). osc_prob_energy_baseline computes the first point in the calling process; if the rest would take under N_JOBS_MIN_PARALLEL_WORK_S (1 s) at that pace, it finishes there too instead of starting workers. A 5-point millisecond scan with n_jobs=2 took 1.43 s and now takes 0.09 s, with the same numbers as n_jobs=1; a scan heavy enough to pay for the workers still uses them. The n_jobs docstring of osc_prob_energy_baseline, which said the argument was forwarded to osc_prob, now says what it does.

• average_spread sets the energy window of average=True across declared discontinuities (issue #134). On that route the window was always the default, ±10%, and a caller’s average_spread was accepted and silently ignored. It is now the half-width of the window as a fraction of the energy, and must lie strictly between 0 and 1. The PhaseAveragingWarning of that route quotes the width and the number of samples actually used, and strategy_info records them. Results at the defaults are unchanged.

• The command-line error for --costhz without a length names every way to give one (issue #136). It said that only --loc-ini/--loc-fin compute the baseline, but --detector-depth and --source-depth do too; the message for --environment earth with neither a direction nor a pair of locations now lists them as well.

• An unrecognized keyword on a wrapper lists the wrapper’s own keywords (issue #152). A misspelled physics keyword (eps_mue, theta12) got a list of engine keywords only; the message now also names the wrapper the caller invoked and its physics and scenario keywords (s12, eps_em, nubar, …), and draws its “did you mean” from both.

Changed

• A short baseline scan at a tight tolerance takes the cumulative scan (issue #125). Under strategy='auto', a single-energy scan of 2 to 7 baselines on a smooth profile at a tolerance tighter than 1e-6 went to the hybrid strategy, since the threshold at which it stands aside, HYBRID_YIELDS_TO_CUMULATIVE_MIN_POINTS = 8, was set at loose tolerances. Below 1e-6 it is now HYBRID_YIELDS_TO_CUMULATIVE_MIN_POINTS_TIGHT = 2. Over 110 scans at rtol = 1e-7 to 1e-12, scored against DOP853, the cumulative scan took 0.06 to 0.18 of the time at the median (at most 1.08, where both sides already ran it), with no silent miss the hybrid did not also make; on realistic profiles the hybrid had mostly declined and the request reached the cumulative scan anyway. Some of these scans now carry a ToleranceNotAchievedWarning, as scans of 8 or more already did, on answers inside the tolerance. Unchanged at 1e-6 and looser, for single points and energy scans, and under strategy='hybrid' or cumulative=False. No result in the paper or its notebooks changes engine. Measurements in docs/dev/measurements/issue125_baseline_scans/.

• cumulative=True over a constant density takes the constant engine (issue #160). A cumulative scan reuses each baseline’s propagator for the next, which a constant Hamiltonian does not need: the constant engine’s one exact exponential per baseline gives the same answer. It ran the cumulative engine, 2 to 3 times slower (0.53 against 0.26 ms for one point, 0.72 against 0.25 ms for 200 baselines, three flavors), and strategy_info now reports 'constant'. A position-dependent potential keeps the cumulative engine, and several energies are still refused.

• diagnostics.rst measures how the true error compares with rtol on smooth profiles (issue #161). The page said an answer that stopped at a tolerance is “typically better” than it, without a number. It now quotes 48 calls on two PREM chords, under three quadrature methods, Magnus orders 2 to 8 and tolerances 1e-4 and 1e-8, against the DOP853 oracle. None was outside the tolerance, the worst was 0.85 times it, and the median 0.014. A reported 1.7x exceedance did not reproduce against that oracle. The script and its output are in docs/dev/measurements/issue161_rtol_gap/.

• Every Hamiltonian builder takes an array for its varying argument (issue #155 §2). The vacuum, *_td, NSI and LIV builders at two to five flavors and hamiltonian_pseudo_dirac_vacuum now accept energy, VCC or l as a number or an array, as the matter builders always have: a number returns one (d, d) matrix, unchanged bit for bit, and an array of shape s returns a stack of shape s + (d, d). Before, an array raised a broadcasting error or, when its length equalled d, returned one matrix mixing the entries without a word. An H_func built from these builders now takes the vectorized path. An array energy is checked entry by entry; scalar calls cost what they did.

• The public arguments are validated once per call, by one set of rules (issue #160). Not every argument in every combination: the checklist audit in docs/dev/measurements/issue160_audit/ lists what is checked, case by case, and the entries above name the cases found open after this one and closed since. Before, most checks lived in the general refinement ladder, so whether a value was refused depended on which engine answered: max_n_slabs=-1 was refused by the ladder and answered by the adiabatic engine. The rules, now applied at the entry of every scenario function, osc_prob_energy_baseline, osc_prob, the public helpers in avgprob, adiabatic, magnus, matter, earth, oscprobstd and plotting, and the command line:

  • real numbers are Python or NumPy reals (so np.float32, np.int64 and 0-d arrays are now accepted), never bool, never complex, and finite; energy and L are checked entry by entry, and L must be at least L0 on every route, the averaged one included;

  • integers are int or np.integer, never bool or a float such as 2.5; flags are True or False only, not truth-tested (average='False' turned averaging on); strings come from their documented set, even where the argument is ignored;

  • rtol and atol are None or above zero; slab and point counts are positive integers; growth factors are above 1 (a factor of 1 never refined); n_jobs is -1 or a positive integer;

  • osc_params, nsi_params and liv_params refuse unknown keys, naming the nearest valid one (a typo’d dcp left dCP unchanged); diagonal NSI couplings, LIV coefficients, sines and phases must be real; Lambda is above zero and n_liv an integer >= 0;

  • t_slab_edges must chain without gap or overlap over the whole path; a user Hamiltonian’s first sample must be a finite, square, Hermitian array; rho_func is sampled at the middle and end of the path as well as at L0;

  • Earth: costhz must lie in [-1, 1] on every route, the geometry helpers included; L may not exceed the chord; locations are (degrees, minutes, seconds) triples in range;

  • the command line refuses non-finite numbers, --precision outside 0 to 17, --dxi13 with --dxicp, and any flag the run would ignore (--rho in vacuum, an NSI coupling without --scenario nsi). Wrong types raise magnus._validate.InputTypeError, which is both a TypeError and a ValueError, so code that catches ValueError keeps working. Messages name the public function the caller called and the argument they passed. The checks never run inside a quadrature or probe loop: arguments still at their defaults are skipped, internal calls skip the helper checks, and single-point and scan timings are unchanged within noise.

• load_nufit_params matches version names without regard to case, like solar-model names (issue #160).

• osc_prob_2nu_matter_std accepts Dm2 = 0, returning the matter-only result instead of dividing by zero (issue #160).

• The 1-3 CP phase is dCP, and the Lorentz-violating one dxiCP, at every flavor count (issue #121). The four- and five-flavor builders called the phase d13 (and dxi13), where the three-flavor ones, the wrappers and load_nufit_params say dCP (and dxiCP), so a loaded parameter set could not be passed straight into a 4ν or 5ν builder. Renamed in mixing_matrix_4x4 and mixing_matrix_5x5, the 4ν and 5ν vacuum builders (all four forms), the 4ν and 5ν LIV builders (both forms), the ten 4ν and 5ν LIV wrappers, the liv_params dictionary of osc_prob_liv, and the command line (--dxicp). The former names still work: d13=, dxi13=, the key 'dxi13' and --dxi13 give the same result, with a FutureWarning; passing both names raises. Positions are unchanged, so positional calls are unaffected. The keyword check costs about 0.1 µs per builder call.

Added

• osc_prob_pseudo_dirac_vacuum, with an exact phase average (issue #165). average=True on a pseudo-Dirac Hamiltonian, through osc_prob_energy_baseline, takes each pair’s phase from the difference of two eigenvalues, which loses a small splitting: against a 50-digit reference it is off by 7e-5 at δm²/max m² = 1e-12 and by 0.25 at 4e-16. The new function takes the pairing (pairs) and the standard parameters, or a mixing matrix and masses, and with average=True forms every phase from the masses and splittings directly, through avgprob.phase_averaged_probabilities_pseudo_dirac; it is within 1.1e-16 of the reference from 1e-4 down to 4e-16. Without average it is the generic route, bit for bit. Nothing that existed changes.

• NSI and LIV oscillograms in compute mode (issue #146 §1). plot_oscillogram, plot_probability_with_profile and plot_biprobability pick the _nsi Earth wrapper when wrapper_kw or osc_params holds an eps_* key, and the _liv one when it holds b*, sxi*, dxi*, Lambda or n_liv. They used to pick the standard wrapper, which refused those keys. Both kinds together, or a flavor count with no such wrapper, are refused with the advice to compute the grid and pass it in.

• nubar= on the probability plots (issue #145 §2). plot_probability_vs_energy, plot_probability_vs_baseline, plot_probability_with_profile and plot_oscillogram take a keyword nubar: in plot-only mode it labels the curves as antineutrino probabilities, which was not possible before; in compute mode it is also passed to the wrapper, so the numbers and the label agree, and a value contradicting wrapper_kw['nubar'] is refused.

• gd.EnergyUnitWarning, for an energy below 1 keV, most likely MeV or GeV left unconverted (issue #141). gd.BaselineUnitWarning now also covers t_breakpoints given in kilometers.

• average_n_samples (issue #134), on every function that takes average: the number of energies the energy-window average samples on a profile with declared discontinuities, 41 by default. The standard error of that average falls as the inverse square root of the number, and each sample costs a full propagation. The other averaging routes sample nothing and refuse the keyword.

• prob is optional on the command line (issue #138). It is the only subcommand, so magnus --flavors 3 ... now runs magnus prob --flavors 3 ..., with the same output; magnus prob ... works as before. magnus --help says so and prints the options of prob after its own, so it lists every flag. The worked examples in cli.rst, installation.rst and the README drop the word.

• hamiltonian_5nu_nsi_td (issue #121), the one builder that existed at two to four flavors but not at five; magnus.hamiltonians now has 43 builders.

• ratio_number_neutrons_to_protons in hamiltonian_4nu_matter_td and hamiltonian_5nu_matter_td (issue #121), a number or a function of position, as in hamiltonian_4nu_matter. Without it, the _td form equalled the constant one only at r = 1. The default, 1.0, is unchanged.

Fixed

• Default axis labels in magnus.plotting no longer print a literal ~ (issue #124). The probability ordinate ('Probability,~', 'Three-neutrino probability,~'), the default profile_ylabel and the baseline and energy abscissae put a LaTeX non-breaking space between text and math, which is a space only when text.usetex is on; under matplotlib’s default mathtext it rendered as a tilde. They now use a plain space, which renders the same under both. A test checks that no label string in the module has a ~ outside math.

• An energy-batched 'gl' scan no longer certifies a sliver of a step onto its slab cap (issue #122). The last refinement is clamped to max_n_slabs, so it can add only a sliver of slabs: on the paper’s five-flavor Listing 1 scan at rtol=5e-13, ladders that had reached 19 926 and 18 346 slabs (floors of 13 and 19) compared them with 20 000 and found agreement to 0.4 and 0.8 of the tolerance; the probability returned was 6.6 times outside it, with no warning. Such an agreement now counts only within the fraction of the tolerance that the step can vouch for, r**p - 1 for a real refinement r at order p (0.015 and 0.41 here), so both scans end in the refinement-caps ToleranceNotAchievedWarning. The capped step that issue #94 keeps certifying (108 -> 115 slabs, agreement 0.14 of the tolerance against 0.28, answer 0.49 of it off) still does. Only the warning changes: across 42 batched scans from 1e-3 to 5e-13, every probability is identical bit for bit; those two scans alone gained the warning. No further level is computed.

• Keywords that were accepted and then ignored now work or raise.

  • default_osc_params_set_name on a two-flavor call (issue #110), where no parameter set applies since sth and Dm2 are required, now raises a ValueError; so does a set name on a scenario function above five flavors, where the vacuum Hamiltonian is given whole. A misspelled keyword on the same call already raised; this one changed nothing in silence.

  • density_matter_is_in_g_per_cm3=True together with density_is_of_number_of_electrons=True (issue #112) now raises a ValueError on all 27 functions that take them. The two give rho different units, and the electron flag used to win: rho=3.0 meant for g cm^-3 was read as 3 eV^3 of electrons, which is vacuum.

  • osc_prob_energy_baseline takes strategy_info (issue #114) and reports the engine on every route, as the wrappers do. It used to reject the keyword as unknown on every route but average=True, which ignored it and, with it, any misspelled keyword; that route now checks its keywords too. Each check is a comparison on a valid call; results are unchanged.

• Arguments other than the energy and the baseline refuse arrays with a message that names them (issue #116). Only energy and L take arrays. An array given for a mixing parameter, an NSI or LIV coefficient, nu_i/nu_f, rho, rho_central, l_scale, costhz or electron_fraction used to fail deep in the call with a NumPy message (“The truth value of an array … is ambiguous”, “setting an array element with a sequence”, “unhashable type”). It now raises a ValueError naming the argument and saying that a scan over it takes one call per value. The checks sit on the path that was already failing, or are a single scalar test, so valid calls are unchanged.

• Energy scans accept a density function that returns a single number, or that takes one position at a time (issue #113). The scenario functions (e.g. osc_prob_matter_std_potential) crashed on an array of energies when given rho_func = lambda l: 3.0*UNIT_G_PER_CM3 (IndexError) or a function using float(l) (TypeError), though both worked energy by energy. The batched engine now broadcasts the first kind, which stays batched, and declines the second to the per-point path. The check reuses the samples the engine already takes, so vectorized densities pay nothing.

• Energy-batched scans with eV-scale sterile splittings no longer lose to one call per energy (issue #111). The batched engine sized its one shared slab grid for the fastest-oscillating energy, so with ~1e4 rad of phase varying as 1/E across a scan it did several times the slab-work of the per-point ladders, and at four or five flavors, where each slab’s exponential dominates, batching came out 2-3x slower than one call per energy. The engine now splits such a scan into groups of energies of similar phase, each on its own grid, choosing the split by a measured cost model (BATCHED_PHASE_GROUPING) from the per-energy phases its seed already computes. Measured on an Earth chord (costhz = -0.5, 40 energies): 4nu with Δm²₄₁ = 1 eV², 0.5-5 GeV, 640 ms -> 220 ms (one call per energy: 370 ms); 5nu with Δm²₄₁ = 1 and Δm²₅₁ = 1.7 eV², 1720 ms -> 570 ms (745 ms); 1-40 GeV, 4nu 320 ms -> 56 ms and 5nu 910 ms -> 130 ms. The split is applied at tolerances of 1e-4 and looser (atol + 0.01*rtol), where the ladders stop on their second level: tighter than that, a group’s coarse first levels could agree by chance (issue #71) where the shared grid had not, and at 1e-8 the shared grid is no more work anyway. Two- and three-flavor scans are never split, and stay bit for bit as they were: there the gain is modest (17-30 % on 3nu scans spanning two to three decades of energy), and a split would give up accuracy those scans have always had beyond the tolerance. A scan the model does not split is computed exactly as before, bit for bit, for at most some 25 us of cost model: the per-energy phases are the singular values the seed always took the maximum of. Warnings stay once per call, and the #71 slab-norm gate applies within each group.

• magnus.magnus_expansion_multislab takes t_breakpoints, with the meaning it has everywhere else: positions where the Hamiltonian is not smooth. The slabs are given here, so the breakpoints are not inserted; they mark which edges are jumps (see Fixed). A breakpoint strictly inside a slab raises UnmarkedDiscontinuityWarning; breakpoints outside the chain are ignored. osc_prob and the cumulative scan pass their own breakpoints down; nothing changes for callers of the wrappers.

• plotting.plot_oscillogram computes the oscillogram when not given one (issue #95). Without a probability array it calls the Earth wrapper (osc_prob_{2,3,4,5}nu_earth) once per zenith angle over the whole energy array, so each call is an energy scan the energy-batched engine answers, with the PREM layer boundaries declared. New keywords: num_flavors, osc_params, the wrappers’ electron-fraction keywords (electron_fraction, electron_fraction_core, _mantle, _crust, _ocean, and ratio_number_neutrons_to_protons; left at None they keep the wrappers’ layered values), wrapper_kw for anything else the wrappers accept, and return_probability. A precomputed probability is drawn as before; giving it together with any of the computing arguments raises rather than ignoring them.

• plotting.plot_probability_with_profile and plotting.plot_biprobability compute through the Earth wrappers too. plot_probability_with_profile takes trajectories (each a costhz, or loc_ini and loc_fin, with its own abscissa and line style) in place of profiles and panels, and makes one wrapper call per trajectory: over its baselines at one energy (x_axis='baseline'), or over its energies at the full chord (x_axis='energy'). Its density panel is then the electron density the wrapper integrates, following the same electron-fraction keywords. plot_biprobability takes configurations (a path, optionally its own energy and parameters, e.g. inverted ordering for one curve) in place of prob_nu and prob_nubar, runs dCP over dcp, and computes markers given as a phase ('dcp'); nu_i and nu_f now also set its default axis labels. Both share plot_oscillogram’s other computing keywords (num_flavors, osc_params, the composition keywords, wrapper_kw, return_probability) and its rule that drawing and computing arguments are not mixed. Existing calls are unchanged.

• strategy_info in vacuum. osc_prob_vacuum and its wrappers (osc_prob_{2,3,4,5}nu_vacuum) accept strategy_info, as every other scenario function does; before, passing it raised ValueError (unrecognized keyword). It reports 'average', 'constant' or 'magnus'. The description of 'declined' now says what it lists: engines that attempted the request and gave up, not every engine that did not apply.

• The phase average, in magnus.avgprob: phase_averaged_probabilities_constant_hamiltonian and phase_averaged_probabilities_adiabatic (issue #64). Every interference term keeps its phase at the central energy and is weighted by exp(-sigma^2 phi'^2/2), phi' = d phi / d ln E, the spread a relative energy spread sigma gives that phase; mixing, crossing amplitudes and the eigenbases at the two ends stay at the central energy. So a phase that runs through many cycles is dropped, as in the L/E -> infinity limit, one that barely moves is kept with its real value, and those between are damped smoothly, where the existing functions keep a pair at zero phase below 2 pi and drop it above (constant Hamiltonian), or drop every pair (smooth profile). On a profile the definition carries into the non-adiabatic windows – an energy offset moves each instantaneous eigenvalue and leaves the eigenvectors – so the answer does not depend on where the windows are drawn: the step on the solar disk of issue #62 came from exactly that. Windows are solved at a few Gauss-Hermite nodes in the offset, the stretches between them are exact phases with exact slopes, and no energy is sampled. Measured against a brute-force average of the same definition on five solar chords from 10 GeV to 10 TeV: within 4.2e-05, where the decohered limit is off by up to 0.14. New constants: AVG_PHASE_SPREAD (0.1, the default sigma), PHASE_AVERAGE_WINDOW_THRESHOLD (0.01, the window search of the profile route), PHASE_AVERAGE_PATCH_ATOL (1e-5, the tolerance of its window solves: at most 4.3e-07 from 1e-7 over 30 solar chords, and up to 18 times faster) and PHASE_SPREAD_SENSITIVITY_THRESHOLD. A pair whose slope is round-off, as a pseudo-Dirac pair’s is over cosmological distances, is kept coherent rather than averaged away on noise (issue #61). The existing functions are unchanged.

• Twelve standard solar models ship with the package, and every Sun entry point – the osc_prob_{2,3,4,5}nu_sun[_nsi|_liv] wrappers and osc_prob_sun – takes density_profile to use one in place of the exponential fit: BP2000, BP04, BS05-OP, BS05-AGS-OP, B16-GS98, B16-AGSS09met, and B23 in six solar compositions, named in any case. The electron density comes from each table’s mass density and hydrogen fraction, interpolated in its logarithm, held flat below the first row and continued along the last interval’s slope past the last. The sterile wrappers take the neutron-to-proton ratio from the same table unless ratio_number_neutrons_to_protons is given; its default is now None, which the exponential fit reads as the 1.0 it always used. stop_at_table_edge=True returns NaN instead, with the new SolarModelRangeWarning, for a baseline that ends past the last row. The new module magnus.solarmodels reads the tables and gives each model’s reference, source and terms of use, and tools/build_solar_model_tables.py rebuilds them from the authors’ files, refusing any whose hash has changed. On the command line, --density-profile takes the model names with --environment sun, and --stop-at-table-edge is new. Defaults are unchanged bit for bit. See the new Standard solar models page of the documentation. Notebook 13 now takes its model by name and compares all twelve on the averaged observable, from a cache (solar_models_cache.json) that continuous integration reads rather than recomputes; notebook 28 draws the solar figures through the Sun wrappers, which moves the 3+1 and 3+2 curves of the averaged-probability figure by up to 3.5e-3, since they now take n_n/n_p from the table.

• average=True on the direct route: osc_prob_energy_baseline, osc_prob_earth and osc_prob_sun now take the keyword that the wrappers have carried since 1.0.0, so a Hamiltonian of your own gets the phase-averaged limit from the same call the wrappers use, by the same three routes – the closed form when the Hamiltonian does not depend on position, adiabatic transport with a Magnus patch at every crossing when it does and the profile is smooth, and an energy-window average when t_breakpoints or t_slab_edges declare discontinuities. The Hamiltonian may be a matrix, a function of the energy, of the position, or of both. average=True refuses return_evolution_operator=True, since the averaged routes form no operator. Defaults are unchanged bit for bit.

• osc_prob refuses average and cumulative by name. Both pass the passthrough guard, because the batching layer declares them, and used to travel down **kwargs to magnus_expansion_multislab, which rejected them with a TypeError naming a function the caller never invoked; the error now says which function takes them.

• Every oscillation-probability function can return the converged evolution operator alongside the probabilities: return_evolution_operator=True makes the call return the pair (P, U), with P exactly what it returns today and U the operator over the same interval, complex and unitary, indexed U[final, initial] so that P == abs(U)**2.T. The keyword lives in the core osc_prob, in osc_prob_energy_baseline and in the generic entry points, and reaches all sixty osc_prob_{N}nu_* wrappers through their **kwargs. It exists for observables built from amplitudes – the mass-state content of what leaves a dense source, the phase-averaged flavor content at a distant detector – which no probability matrix can give back, since the phases are gone. With the keyword set, the refinement ladder compares the operator itself between levels, at the same rtol and atol, so the returned operator is converged in its phases and not only in its moduli; the specialized engines (hybrid, interaction picture, scan) stand aside for the call, since only the general ladder forms the operator, and a baseline scan takes the per-point path instead of the cumulative traversal. average=True and strategy='hybrid' are refused together with it, with an error naming the entry point, rather than ignored. Every default result is unchanged bit for bit: with the keyword off, the only difference is one boolean test at the exit points of the core.

• Either end of an Earth trajectory can now be underground. Every Earth entry point takes source_depth and detector_depth, and earth.distance_traveled_inside_earth, earth.earth_radial_distance_from_depth and earth.prem_layer_edges_along_chord take them too. The zenith angle is measured at the detector, which is what it already meant when the detector was on the surface, so a buried detector also sees downward-going neutrinos (costhz > 0) through its overburden — a trajectory the surface geometry gives no path for at all. Naming detector_depth fixes where the trajectory ends, so L is then computed rather than given, and passing both raises. The two named locations still describe a surface-to-surface chord, so combining them with a depth raises as well. Both defaults are zero and every default result is unchanged bit for bit: the three geometry functions return through the expressions they have always used, rather than through the generalized ones, which agree on every zenith angle tested but are not the same expression.

• density_matter_ocean replaces the density of PREM’s outermost shell, wherever electron_fraction_ocean already replaced its composition. That shell is 3 km of global-average ocean at 1.020 g/cm³; a detector under continental rock sits under about 2.6 instead, and one under Antarctic ice under about 0.92. It is a correction worth making for a trajectory close to horizontal, which can spend its whole length inside that shell.

• Every NuFit release now has a named oscillation-parameter set, so default_osc_params_set_name reaches all eighteen of them rather than only 6.0 and 6.1. From 4.0 onward a release splits its fits by whether Super-Kamiokande atmospheric data is included, and both halves are named: OSC_PARAMS_NU_FIT_<version>_SK_<ordering>, the spelling the 6.0 and 6.1 entries already used, and OSC_PARAMS_NU_FIT_<version>_NOSK_<ordering>. An earlier release has no such split and is OSC_PARAMS_NU_FIT_<version>_<ordering>; its secondary categories stay reachable through load_nufit_params alone, since NO already means normal ordering here and a name carrying huber_fluxes_no_rsbl would read as two orderings at once. Fifty-three names in all. They are generated from NUFIT_GLOBAL_FITS, so a future release needs no second table. The four existing entries are untouched: the 6.0 dictionaries are built from module constants whose values differ from the loader’s in the last bit, and a caller pinned to 6.0 keeps the same bits.

Changed

• average=True states the neutrino’s initial state, the flavor state by default (issue #73). A new keyword, average_initial_state ('flavor' or 'decohered'), on every function that takes average. The initial state used to depend on the route: the constant route started in flavor, the smooth route decohered. The same medium therefore gave two answers: P_mumu = 0.3626 at 1 GeV, 300 km and 3 g/cm^3 on the constant route, and 0.4528 as a flat exponential profile. Now both give 0.3626 by default, and 0.4528 with average_initial_state='decohered', on both routes.

  • What changes: on a smooth profile, a neutrino produced in the medium keeps the interference between eigenstates present at production, damped by the spread of its phase from there on (the paper’s Eq. (flavor_vs_decohered)). With no non-adiabatic window this is a closed form, Eq. (phase_average_nocross). With windows, the coherence is carried from the start of the path into the first window. avgprob.phase_averaged_probabilities_adiabatic takes the initial density matrix as rho0, and its default is now the flavor state.

  • What does not: a neutrino produced in the solar core has every such term damped from production on, and the solar MSW curves are unchanged, bit for bit (tested on BS05-AGS-OP at 40 energies). So are the paper’s geoneutrino averages and its supernova curve, at 59 of 60 energies. Neutrinos that enter a medium from outside, as in solar tomography, need average_initial_state='decohered' to reproduce the previous numbers; with it they do, bit for bit. Notebook 28’s solar-disk chords now pass it.

  • Declared discontinuities: that route averages over an energy window from the flavor state, as before; 'decohered' raises ValueError there.

  • Speed: before transporting a flavor start, the phase slopes are estimated on 17 points, about 80 us an energy. Where every interference term would be damped anyway, the decohered start’s value is returned. With the phase average on, the smooth route also stops computing the coherence report it never reads, about 0.2 ms an energy. The net effect: the 40-energy solar curves are 2-9 % faster, and nothing measured is slower beyond noise.

  • Memory: the phase average now bounds its interference terms by memory (4 GiB) as well as by number, and falls back to the limit with PhaseAveragingWarning where the bound is reached. At five flavors on a three-window profile, the decohered start used to exhaust 12 GB before reaching the count bound, on main as well.

• The energy-batched engine seeds its starting slab count for 'trapezoid' and 'simpson' too, when the phase-based seed is at least QUADRATURE_SEED_MIN_SLABS (4); below it they start from min_n_slabs as before, and 'gl' is unchanged. The threshold is measured: over 212 scans (2 to 40 energies, smooth and breakpoint profiles, both rules, rtol = atol = 1e-3 and 1e-4), a seed of 2 or 3 could cost up to 2.14x the work, and a seed of 4 or more never cost any (at 4 the worst case did 0.87x the work; at 5 and above, 0.36x or less). On a core-crossing chord at 50 energies, 'simpson' at 1e-4 went from 1016 ms to 40 ms. The per-point ladder keeps its 'gl'-only seed.

• Notebook corrections from a read of all twenty-nine notebooks (notebook 28 left for its own rerun). Wrong numbers: 05’s vacuum biprobability used 810 km and 2 MeV under a T2K label (now 295 km, 0.6 GeV); 01 compared a neutrino vacuum curve with antineutrino matter curves, and in the energy plot at a different baseline; 16, 17 and 19 built PREM with a uniform Y_e = 0.5 and now use the layered composition the Earth wrappers use (16: the mean-density substitution is wrong by up to 0.59, was 0.51; 17: the core-chord ordering split is 0.43 / 0.42, was 0.49 / 0.42). Text that contradicted the library: 01’s two- and three-flavor Hamiltonian equations, its claim that the three-flavor probability matrix is symmetric, a function name that does not exist, and osc_prob taking an array of baselines; 11’s order ceiling (10, not 6); 24’s Gauss-Legendre ceiling (8, not 6); 25’s four-flavor exponential (a Jacobi eigensolver since 1.0.9); 20’s warning count (fifteen) and false-alarm attributions; 21’s refinement fix (PR #94 alongside #35); 29’s statement on the regime between the averaging limits. Text that contradicted its own output: 12’s comparison of the three strategies (re-run and rewritten), 24, 25 and 27 (whose ladder table counted MagnusConvergenceWarning as running out of room). Also the broken Mag$\nu$s in the 10 and 12 introductions, stray return values printed after plots, 14’s section numbering, 09’s listing of the _liv wrappers, 19’s “10%%”, and a note on the expected tolerance warnings in 02’s and 03’s solar sections.

• Notebooks 04, 05, 06 and 07 compute their Earth probabilities through the Earth wrappers (issue #95 for notebook 06). They built the PREM Hamiltonian by hand, without declaring the layer boundaries, on fixed grids (notebook 06 also with a uniform electron fraction of 0.5). The wrappers declare the boundaries, use the layered electron fraction and refine to their default tolerance. One exception is stated where it happens: notebook 07’s 3+1 energy scan reaches the standard max_n_slabs at its lowest energies, where the sterile phase is of order 1e6 radians, so it warns as expected; a new cell computes one of those probabilities with max_n_slabs=10**6 as well, which removes the warning at a longer runtime. Against a reference at rtol=1e-9, the largest error on the sampled points fell from 3.2e-2 to 8e-12 (07, 3+1 to the South Pole), from 1.2e-3 to 5e-6 (05, 2 GeV to the South Pole) and from 3e-4 to under 1e-6 (04). Notebooks 04 and 05 now call the computing forms of plot_probability_with_profile and plot_biprobability; notebook 07’s 3+1-versus-energy figure labels its axis in MeV, the unit its energies were always in, and its aliasing comment gives the sterile phase as ~1e6 radians (it said ~1e9). Notebook 06’s two-flavor nu_mu-nu_tau oscillogram is now computed in vacuum along each chord: neither flavor feels the charged-current potential, so the matter term drops out, and the two-flavor Earth wrapper (which applies it to its first flavor, nu_e) had given this system a matter effect it does not have. Its crude-settings demonstration moved to the 1-2 sector, with its error stated (up to 2e-2 at one slab per layer and first order, against 1e-5 at the default tolerance). Notebook 05’s best-fit markers are labeled NuFIT 6.1, the parameters they use (they said NuFit 6.0), and notebook 01’s text says the same of its parameters. Notebook 05’s check against the closed-form vacuum formula now draws its own markers: it drew the Magnus markers of the cell before, and its unused inverted-ordering markers had been computed with the normal-ordering mass splittings and best-fit phase.

• Python 3.13 declared supported. The pyproject.toml classifiers now list 3.13, which CI has tested and passed on main; they listed only 3.10-3.12 before, so PyPI did not show 3.13 as supported.

• strategy='auto' hands a moderate phase to the Magnus ladder (issues #70 and #120). On a smooth profile it used to run the hybrid strategy at every tolerance, and the hybrid’s cost is its window search, which does not follow the tolerance: the four scans of the paper’s Fig. 1 took 8 s at the default of 1e-3. Now the ladder answers instead when min(rtol, atol) is at least AUTO_LADDER_MIN_TOLERANCE (1e-6), the estimated accumulated phase (the integral of the spread of H’s eigenvalues) is at most AUTO_LADDER_MAX_PHASE (1e4 rad), and the ladder can start within a quarter of its slab cap (AUTO_LADDER_MAX_FLOOR_FRACTION). It runs at a tenth of the tolerance (AUTO_LADDER_TOLERANCE_MARGIN), without the interaction picture, starting on slabs over which the Magnus series is guaranteed to converge. Over 20 smooth workloads with phases from 5 to 1.2e4 rad it was 2 to 60 times faster on a single point and 12 to 500 times faster per point of a 40-energy scan, within the tolerance on every one; the Fig. 1 scans now take 40 ms of computation. strategy_info reports the handoff as the hybrid declining, with the reason 'auto prefers the ladder', and the hybrid’s test for an undeclared density jump still runs, with its reason and UnmarkedDiscontinuityWarning. An energy scan at one baseline that the energy-batched scan will take goes to it whatever its phase (issue #84): the phase limit prices the ladder point by point, and applied to a scan it sent one three-flavor, two-resonance scan of 100 energies to the hybrid, 62 s against 0.4 s. Below 1e-6 (issue #120) the route stays open on integration_method='gl' at a single baseline, for a phase within a limit that shrinks as (tol/1e-6)**(1/p) with the requested order p, capped at AUTO_LADDER_TIGHT_MAX_PHASE (2 000 rad) and applied to energy scans as well; the ladder then runs at the tolerance itself. This is what lets the paper’s Listing 1 (rtol=1e-12, atol=1e-14, magnus_exp_order=8) drop strategy='magnus': its four curves go to the ladder, which answers them in 0.015 to 0.13 of the hybrid’s time, 2.3e-14 from DOP853 at the worst point. Over 139 workloads at 1e-7, 1e-9 and 1e-12 the ladder at order 8 missed no tolerance without a warning and never warned where the hybrid had certified. At order 4, four workloads at 1e-12 came back 1.01 to 1.15 times outside it; the hybrid misses four there too, by up to 1.3 times. Every solar path, baseline scans and other quadratures below 1e-6, an invalid order (same error), strategy='hybrid' and strategy='magnus' are unchanged, and no notebook’s engine changes.

• average=True returns the phase average (issue #64), with the spread set by a new keyword, average_spread (default 0.1), on every entry point that takes average. Before, a constant Hamiltonian kept each pair of levels at zero phase below 2 pi and dropped it above, and a smooth profile dropped every interference term, including the readout at the end of the path: so a phase of 1 rad was set to zero, one of 20 rad was dropped although 14 per cent of its interference survives a 10 per cent spread, and the solar disk stepped wherever the window search changed its windows (issue #62). Each point is still computed the old way first, and the old value is returned bit for bit wherever the phase average agrees with it to 1e-4 – on a profile it is not even recomputed without a non-adiabatic window, since adiabatic transport of a decohered start carries no interference. Measured over every averaged call in the notebooks, tests and documentation, about 5 900 points: every solar MSW curve, notebooks 13, 24, 25 and 27, and paper Figures 14-17 are unchanged; 1 557 points move, 1 538 of them pixels of paper Figure 5f from 10 GeV up, the rest short baselines and profiles with windows at GeV energies, each a case where the old value was not the average. PhaseAveragingWarning now says that the result depends on the spread (it changes by more than 1e-3 per e-fold of it) rather than that no averaged expression applies. A Hamiltonian that does not depend on energy, passed as a matrix or as a function of position alone, has no spread to average over and keeps the old behavior; so does the energy-window route of a profile with declared discontinuities. strategy_info records the spread, how many points were recomputed, and the largest sensitivity. Where the result does not change it does not cost more either: on a constant Hamiltonian the points whose pairs have all decohered are now formed in one batch, bit for bit, so an astrophysical vacuum average over 1000 energies takes 11 ms against 26 ms before, check included.

• A converged average=True no longer warns about the coarse levels of its windows (issue #66). The phase average solves each non-adiabatic window on a ladder of slab counts, and its first, coarsest levels raised MagnusConvergenceWarning although the level returned was fine: the chord of the paper’s solar-tomography listing printed it four times. The window ladders now check the slab norm once, on the level they return; the hybrid strategy and the decohered route keep the check on every level. Over 55 solar chords, firings fell from 6 to 3, none of those dropped preceded an error above 1e-4, and the probabilities are unchanged bit for bit.

• rtol and atol set the tolerance of the phase average on a smooth profile (issue #65). They were accepted there and dropped without a word: the window patches and the stretch phases converged to a fixed 1e-5, and PHASE_AVERAGE_PATCH_ATOL, bound as a default argument at import, did nothing when changed. Both now converge to the tighter of rtol and atol, 1e-3 by default, and the decohered limit is returned only where the phase average agrees with it within that tolerance, or within 1e-4 if that is tighter. Over twenty chords through the solar core, 30 GeV to 3 TeV, the default moves the probability by at most 4.6e-06 from its value at 1e-5 and is 2.4 times faster at the median (0.9 to 6.3). phase_averaged_probabilities_adiabatic takes patch_atol and phase_tol; their defaults, PHASE_AVERAGE_PATCH_ATOL and the new PHASE_AVERAGE_PHASE_TOL (1e-5 each), are read at each call.

• The adiabatic machinery evaluates the Hamiltonian in batches. The search for non-adiabatic windows took the finite-difference derivative at every probe point, and bisected every gap extremum, one Python call per position: about 1300 Hamiltonian calls per energy on an averaged solar call, most of its time. The probe grid, the bisections (all of them advanced together, sixty batched steps), the adiabaticity at every candidate, the growth of every window, the parallel transport of adiabatic_propagator, the accumulated phases and oscillation_sampling now evaluate the Hamiltonian over an array of positions and decompose it in one batched call. Every reduction is still formed per position as before, so the results are unchanged: 15 044 calls to the adiabatic and averaging functions, over the four solar curves, disk chords with windows, hybrid-strategy calls, notebooks 10, 12, 13 and 23 and the adiabatic and averaging tests, return identical bits. The one difference found was a single derivative entry, 1 ulp, on a chord profile whose scalar and array evaluations disagree in the last place (NumPy squares an array but calls pow on a scalar); nothing downstream moved. A Hamiltonian that only takes one position at a time is evaluated position by position, as before. The four averaged solar curves of the paper’s solar figure (three flavors, with NSI, 3+1 and 3+2, ninety energies each), timed interleaved against the previous code on an idle machine: 1.6, 1.6, 3.1 and 4.9 s before, 0.56, 0.57, 0.88 and 1.2 s after, 3.5 times faster together. The hybrid strategy shares these functions and gains with them (issue #64).

Fixed

• The 3+1 and 3+2 mixing-matrix docstrings now describe the matrices correctly. mixing_matrix_5x5 said it carried three active-sterile CP phases (δ14, δ15, δ35), but it carries four (δ14, δ15, δ24, δ35). It now also explains why it has nine angles and five phases: the missing pair belongs to the rotation between the two sterile flavors, which changes no probability among the active flavors. Both docstrings credited arXiv:1103.4570 to Kopp, Machado, Maltoni & Schwetz; its authors are Kopp, Maltoni & Schwetz. The companion reference arXiv:1105.3911 now names its authors, Barry, Rodejohann & Zhang. The matrices themselves are unchanged.

• The hybrid strategy no longer certifies a wrong answer across a sharp step between its probe points (issue #100). A density step narrower than the probe spacing, away from any resonance, is neither a gap extremum nor a probe point, so no adiabaticity check ever landed on it: every refinement agreed on adiabatic transport across it, and strategy='auto' returned the hybrid’s answer as certified, with no warning. On the issue’s profile (a shock of width 0.8 on a path of 3e5, in units of the inverse vacuum splitting) that was P_ee = 0.153 against a correct 0.492. Before certifying, hybrid_propagator now compares, on every probe interval, the change of H across it with the derivatives measured at its ends, which catches a step narrower than the probe spacing wherever it falls between two probe points (one within a few widths of a probe point is seen by that point’s own derivative). At the steepest point of each interval flagged it evaluates the adiabaticity parameter exactly, and if that could move the answer by the tolerance it repeats the refinement with the point examined from the first level. A restarted pass looks again, on the finest probe grid too, past the points it already examines, for up to 8 restarts, so a profile with more steps than one look examines (16) is covered as well: a staircase of 32 small steps, each harmless alone, went from certified 1.8e-02 off to 1.1e-04. The issue’s case now certifies 0.4921 against 0.4922.

  • Nothing else moves. Over the 53 cases of the physical-profile population (solar tables, BS05, supernova shocks and turbulence, Earth chords with an undeclared crust, and the issue’s profiles), 50 are bit-for-bit what they were, probabilities and certification alike. The three that change were certified wrong by 0.34 to 0.49, and now certify within 1.2e-4 of solve_ivp.

  • Nor does the speed. The check reads arrays the probe sweep already computed: 25 us at 400 probe points and 0.31 ms at 6400, only on calls about to certify. Only a call that would otherwise have certified a wrong answer pays for a restart. Over 37 timed cases, from 7 ms tabulated profiles to 5 s solar-model scans, main and this change agree within run-to-run noise (-1.1 % summed).

  • It still cannot see a bump that rises and falls between two probe points (see docs/source/adiabatic_strategy.rst); a step, however narrow, is no longer in that class.

• The energy-batched 'gl' scan no longer accepts an energy whose own slabs are still wide (part of #71). On smooth profiles two coarse levels could agree by chance while each slab still spanned several radians of phase, and the scan returned that agreement silently: over the 477 scans of the #71 pool, 27 energies up to 19 times outside the tolerance; with antineutrino and magnus_exp_order 2 and 6 variants added (191 'gl' scans), 72, up to 102 times outside it. The scan now refuses the agreement of every energy whose largest slab norm is BATCHED_GL_MAX_SLAB_NORM (2 pi) or more, while it can still add slabs, and that energy goes one level further. The 27 go to 0 and the 72 to 2, with none added, for 0.9 % and 3.1 % more total work. The norms are the eigenvalues _expm_stack already computes (magnus._row_slab_norms), and MagnusConvergenceWarning is unchanged. Every scan the gate leaves alone is bit-for-bit what it was, and so is every other path, at no measurable cost. Where the phase is concentrated, the gate can refuse correct agreements for several levels: 9 of the 191 scans that had been within tolerance do more work, the worst a 40-energy antineutrino scan at rtol = atol = 1e-3 (3.2x to 3.6x the time, 6 ms to 20-23 ms). The alternatives measured to avoid that cost all left more misses or cost more (see the constant’s docstring).

• The energy-batched scan now warns when it accepts an energy at the slab cap (part of #71). Once 'trapezoid' or 'simpson' reach max_n_slabs (2000 by default), each further level refines only the points per slab; two such levels agreeing verifies the quadrature inside each slab but not the slab count, and the scan accepted them silently. On the multi-resonance test profile at rtol = atol = 1e-6 that answer was 2.5e-4 off (250 times the tolerance) while successive levels agreed to 7e-8. Such an acceptance now raises ToleranceNotAchievedWarning, naming max_n_slabs and integration_method='gl'. Nothing else changes: the same levels, work and probabilities. Over the 477 scans of the #71 measurement pool it catches 62 of the 132 energies that were outside tolerance without a warning; it also fires on 297 energies that happened to be within tolerance, all in the three workloads that need more than 2000 slabs, where their accuracy could not be verified either.

• 'trapezoid' and 'simpson' now keep their order at declared breakpoints (part of #71). Both sample each slab at its two ends, and a declared breakpoint is a slab edge; at a density jump the profile returns one side’s value there, so the slab on the other side integrated a sample from the wrong layer, an error falling only as 1/n_tpts_per_slab whatever the rule. That sample is now taken just inside its own slab (by min(max(1e-8 w, 8 ulp), w/2), w the slab width), in the per-point ladder, the cumulative scan and the energy-batched engine. 'gl' samples no endpoints and is bit-identical. Measured on PREM chords and the Earth wrappers at rtol = atol = 1e-3 and 1e-4: silent misses of the tolerance fell from 122 to 19 in the energy-batched engine (those left were accepted at a slab norm of pi or more, a separate stopping-rule question) and from 19 to 0 in the per-point ladder, with the remaining quadrature error falling at the rule’s own order; the cumulative scan at its default 100 points per slab went from 4.9e-5 to 1.6e-12 ('simpson') and from 1.8e-5 to 1.5e-9 ('trapezoid') on a core-crossing chord.

• The energy-batched engine could stop on two nearly identical grids and return, without a warning, an answer outside the requested tolerance, on any energy scan with declared breakpoints: Earth chords through the Earth wrappers or t_breakpoints (issue #71). The breakpoints are inserted into every refinement level’s grid, so at small slab counts a nominal x1.5 step barely refines it: on the PREM chord at cos(theta_z) = -0.9, 4 -> 6 slabs is 20 -> 22 points. At two flavors and rtol = atol = 1e-3, two of twelve energies stopped on that comparison and came back 2.5e-3 off, and the three-flavor Earth wrapper at 1e-4 came back 2.8e-4 off. The engine now applies the per-point ladder’s bound, MIN_EFFECTIVE_REFINEMENT, on grids with breakpoints, and grows the real grid instead of the nominal slab count, so it does not spend levels on steps the bound would refuse. Every case measured is now within tolerance. Scans are faster at two flavors (about 0.7 of the time) and at three flavors from rtol = 1e-4 on (0.7 to 0.9); at three flavors and 1e-3 they take 1.1 to 1.2 times as long, since the old ladder stopped there on a false agreement. Without breakpoints, and at four and five flavors on the PREM chord, results are unchanged, bit for bit.

• The hybrid strategy could certify, without a warning, an answer outside a tolerance tighter than about atol + rtol = 1e-6. hybrid_propagator checked its result against the requested rtol/atol, but every Magnus patch inside a non-adiabatic window converged only to a fixed 1e-7, and the agreement test cannot see a patch’s error: a window that does not move between refinement levels holds the same patch in both. The paper’s Listing 1 at rtol=1e-12 (3nu, 2 MeV) came back 2.0e-11 off DOP853, and a 2nu solar chord (0.3 R_sun, 10 MeV) 4.8e-09 off at every tolerance, ten times outside rtol=1e-9, atol=1e-11, certified each time. Every patch now converges to min(1e-7, (atol + rtol)/10), so at atol + rtol >= 1e-6, the default included, nothing changes, bit for bit. Over 140 cases at rtol 1e-7, 1e-9 and 1e-12 checked against DOP853, 78 were silent misses under strategy='hybrid' and 'auto', and none are now. A patch that cannot reach the tolerance within its slab cap leaves the result uncertified, with HybridCertificationWarning, as before: that solar chord is now 3.0e-10 off, and warns at rtol 1e-9 and 1e-12.

• MAGNUS_PAPER_CACHE_ONLY now forbids notebook 28 from recomputing anything, as both READMEs said it did (issue #63). Only the scan and timing sections honored it; cached(), the helper behind about thirty others, printed “configuration moved, recomputing” and recomputed on the runner, so a section whose configuration moved passed continuous integration instead of failing and naming itself. It now stops the build the way those two do. The committed cache holds every section under its current key: rebuilt with the variable set, all 84 are read back and none recomputed.

• average=True missed any feature narrower than the grid it searches for non-adiabatic windows on, and returned the fully adiabatic answer without a warning (issue #60). The adiabatic averaging engine looks for windows once, on 200 probes, and never refines; on a supernova shock ray 70,000 km long the probes are 350 km apart, so fronts 0.07 to 70 km wide were never examined, and the call returned 0.04 where the averaged probability is 0.37 to 0.59. The profile is now checked first for features that sharp and able to move probability between levels – an instantaneous change across them would move more than the new avgprob.SUDDEN_TRANSFER_THRESHOLD, the default tolerance. Where there is one, the windows come from the hybrid strategy’s refinement instead, which certifies them; where no refinement resolves the feature, UnmarkedDiscontinuityWarning says so and names the cure, t_breakpoints. Of 24 fronts 0.07 to 2000 km wide on that ray, 16 are now within 0.01 of a decohered reference and 8 warn; 16 were silently wrong before. Nothing else moves: every Sun wrapper, the BS05 solar model and the tabulated profiles return the same result, bit for bit, and of the roughly 5,000 averaged calls in the notebooks only 16 pixels of paper Figure 5f are escalated, three of them changing, each to within 0.001 of a reference. The report of averaged_probabilities_adiabatic gains escalated, resolved and certified. The issue’s own reference values, 0.84 and 0.18, were a position average at one energy, which keeps the interference between the two fronts; the averaged probability is the decohered one, and against it the answers with the fronts declared are right, not wrong as the issue said.

• The hybrid engine carried a decoupled state across an exact level crossing onto the other level, and certified the result (issue #59). At their defaults the 4nu Sun wrappers decouple the sterile state, whose matter term then crosses the lowest active level inside the Sun; for neutrinos between 1 and 300 MeV, osc_prob_4nu_sun, _nsi and _liv came back certified and wrong at 17 of 25 energies, the sterile row landing on an active level (P_ss = 0 where it is 1) and P_ee off by up to 0.53. Two levels that do not couple have an adiabaticity parameter of 0 over a gap of round-off, and the test for a degenerate pair was gap > 0, which caught an exact crossing only when round-off happened to give a gap of exactly 0.0. A gap eigh cannot resolve now counts as degenerate, so every such crossing gets a window; the bound is the new adiabatic.DEGENERACY_ULPS, measured at 303 exact crossings against 2565 ordinary candidates, with eleven orders of magnitude between them. Every call without a degenerate pair is unchanged by this part, bit for bit.

• The hybrid engine’s refinement certified on agreement alone as soon as any window had opened, so a window at one resonance vouched for the stretch of path it did not cover: a resonance just below the adiabaticity threshold elsewhere was left to adiabatic transport, which agrees with itself whether or not it is right. The requirement that the adiabaticity parameter fit the tolerance, which used to apply only when no window opened, now applies to whatever no window covers. Found through the fix above, where the window at the sterile crossing let osc_prob_4nu_sun at 237 MeV certify 1.6e-03 out after two iterations while the 3nu call on the same physics refines five times and is right to 8e-05. Over 300 Sun calls of two to five flavors this moved eight results besides the 4nu defaults: seven became 55 to 7500 times more accurate (the largest gain, a 3nu call at 187 MeV, from 1.3e-04 to 1.8e-08), and one the hybrid can no longer certify is answered by the general ladder within the requested tolerance. find_nonadiabatic_windows and hybrid_propagator report the new quantity as gamma_unpatched in info.

• The scenario functions warn, past MAGNUS_MAX_PREDEFINED_NUM_FLAVORS, that they will use the vacuum Hamiltonian passed in h_vac_energy_indep, and that path did not work. unpack_oscillation_params_from_dict fell off the end of its own branch and returned None, which validate_input_battery then iterated, raising TypeError; with validate_input=False the parameter-filling step ran on names that the 2-to-5 unpacking had never assigned, raising UnboundLocalError. The helper now returns an empty array, and the four parameter-filling sites are bounded above as well as below, so osc_prob_vacuum, osc_prob_matter_std_potential, osc_prob_matter_nsi and osc_prob_liv accept any flavor count when handed a Hamiltonian. Verified at six and eight flavors, in vacuum and in matter; the two- through five-flavor results are unchanged, bit for bit.

• earth.dms_to_decimal added the minutes and seconds to a negative degree value instead of counting them in its direction, so every West or South coordinate with non-zero minutes landed up to a degree too close to the meridian or the equator. The function sits under chord_length_inside_earth, costhz_between_points_on_surface, every osc_prob_*_earth call made with loc_ini/loc_fin, and the CLI’s --loc-ini/--loc-fin: the chord from Fermilab to Homestake came out 1207 km where the true chord, and DUNE’s baseline, is 1285 km, with the trajectory’s direction wrong by the same amount. Sites east of Greenwich and north of the equator, and the poles, were unaffected. The sign is now read from the first non-zero part and the minutes and seconds are taken as magnitudes, so (-88, 15, 26) and (-88, -15, -26) both read as -88.257 degrees. Regression tests added.

• Documentation brought in line with the code, from a read of every docstring and of the README and the documentation pages against the code they describe. The notable corrections: a recipe that raised (d13= for dCP=) and abs(U)**2.T, which does not parse; the dispatch order under strategy='auto', which is not silent when the hybrid declines an undeclared jump, and does not try the interaction picture first; average=True described as the decohered limit where it returns a phase average; the runtime dependencies (five, not three: matplotlib and numba were missing, as in src/requirements.txt); the location format of the Earth wrappers; the matrix-exponential backends, unitary to round-off rather than by construction; stale counts, Returns/Raises sections and cross-references; two docstrings that contradicted the provenance notes beside them; versionadded:: 1.0.5, a release that does not exist (1.0.11); and CITATION.cff, which still named 1.0.0.

[1.1.0] - 2026-09-06

Changed

• The sterile block of the matter projector now follows the Earth’s layered composition. The matter Hamiltonian factorizes as V_CC(l) times a projector, and the density has always resolved the electron fraction per PREM layer – but the projector was one constant matrix taking a single scalar r = n_n/n_p for the whole trajectory. The two therefore described different media on any chord whose composition varies, which is every chord that reaches the core.

  Only the sterile states are affected: the active flavours all feel the same neutral-current potential, so it drops out as a phase, while a sterile state feels neither current and is left carrying -V_NC = (r/2) V_CC. Three flavors are therefore untouched, exactly – the projector’s sterile block is empty there and the constant matrix is still returned.

  The cost of getting this wrong was larger than the shipped warning said. Off resonance it is the ~2e-02 the warning quotes, but near the sterile matter resonance on a core-crossing chord it reaches 0.4 in probability, flat under refinement, invariant in the mass splitting, surviving small mixing angles, and present at 3+2 as well. No choice of scalar fixes it: the best value obtainable by scanning still leaves 7e-03 in that band.

  It costs nothing to fix. Every Earth entry point already declares the PREM layer boundaries as slab edges, automatically, so each segment is homogeneous in composition and a per-segment projector is exact on the grid that exists – no regridding and no extra Hamiltonian evaluations. Measured at parity with the constant-projector path.

  This changes what an Earth call returns by default at four and five flavors. ratio_number_neutrons_to_protons defaults to None on the twelve Earth wrappers, resolving to the layered composition; passing a scalar explicitly reproduces the old behavior exactly. No signature changed: the parameter accepts a callable as well, and the general (non-Earth) entry points keep their scalar default, having no layered composition to resolve.

  SterileMatterCompositionWarning is retargeted rather than removed. The default and a caller-supplied callable are silent; a scalar over a layered profile warns unconditionally, the previous 2% threshold having been measured passing chords whose error equalled the figure the warning quoted.

  notebooks/sterile_projector_check.py reproduces the defect and the fix in one command, three arms differing only in the projector.

[1.0.13] - 2026-09-06

Fixed

• The HiddenFeatureWarning told callers to pass t_breakpoints and then printed ones that did not work. The suggestion was a pair-scale bracket, ten to a hundred times wider than the feature it was meant to straddle: a single point stayed at 3.0e-02, and a 60-point scan moved to 5.8e-02 with no warning at all, because passing t_breakpoints switches the scan the detector runs in off. A caller who did what the message said therefore landed back where the warning exists to rescue them from, with a worse error and nothing on screen.

  The suggested edges are now localized: the flagged interval is re-sampled and the sub-interval carrying most of the variation is bracketed instead. Following the message takes the same cases from 3.0e-02 to 1.0e-04 and 1.0e-02 to 9.9e-05, at a point and over a 60-point scan alike. The localization runs only where the warning fires, so the detector’s false-positive record is untouched by construction.

  The message also now says what to expect afterwards: breakpoints route the call to the general slab ladder, so accuracy lands near the requested tolerance rather than below it.

  Detection itself is unchanged. The undetectable remainder is structural – a feature narrower than the reference grid cannot be found by refining a grid – and no attempt was made to change that.

[1.0.12] - 2026-09-05

Changed

• The baseline scan folds its running product in a compiled kernel rather than a Python loop. :func:osc_prob_energy_baseline answers every requested baseline from one traversal of the profile, because each answer is a prefix of the next; the loop that walked it snapshotted the running product at each requested distance, one Python iteration per slab. Marginal cost per slab falls from about 1.2 us to 0.011 us at two flavours.

  End to end the gain grows with the grid, since the fold’s share of the scan does: about 2.7x at 512 accumulation steps and 7.5x at 32 768 at two flavours, 1.5x and 3.9x at three. It exceeds the survey’s own 2-4x estimate for a reason worth recording – the kernels of 1.0.6 through 1.0.11 had made building the operators so cheap that the Python fold had become the dominant cost of this scan.

  This needed its own kernel rather than reusing 1.0.11’s. That one folds acc <- acc @ U[k] with k descending; this one needs U[k] @ acc with k ascending, which is the mirror image – a different product, not a different parenthesization – and it must snapshot mid-fold at requested baselines, which a final-answer kernel cannot do.

  Not bit-identical, for the same reason as 1.0.11: NumPy routes these products through MKL’s zgemm, which uses fused multiply-add, and a kernel compiled with fastmath off cannot reproduce that ordering. Worst shift 1.4e-14, which is twenty times smaller than 1.0.11’s. Unlike that release, neither side is systematically closer to exact – against 40-digit mpmath the new fold wins at two flavours and the old at three – so this is rounding exchanged rather than accuracy gained. A numba-less install is exactly bit-identical to the old loop. The scan-against-per-point invariant was checked on identical slab edges at all four flavour counts, 4.9e-15 to 2.2e-14, against the suite’s 1e-12 bar.

[1.0.11] - 2026-09-05

Changed

• The interaction-picture engine folds its slab operators in a compiled accumulator rather than a Python loop. That loop ran one iteration per slab, each doing a single small matrix product and returning to Python, and measured 56% of the engine’s pass; the fold itself falls 18.3x, the engine 2.28x. What a caller sees depends strongly on how many slabs are folded – about 2.3x at 32 768 slabs, and nothing measurable at 512 – so the figure is a property of the request, not of the engine.

  This is the fourth loop of its kind to be compiled, after the separable scan’s in 1.0.6 and the Gauss-Legendre expressions in 1.0.8, and it reaches a narrower set of calls than any of them: the engine serves only two-flavour requests on a profile tagged as exponential, with no supplied slab edges.

  Not bit-identical, and established as unachievable before it was written rather than discovered afterwards: NumPy routes these stacked 2x2 products through MKL’s zgemm, which uses fused multiply-add, and a kernel compiled with fastmath off cannot reproduce that ordering. Worst engine-output shift 2.8e-13 over 34 configurations spanning slab counts 8 to 32 768 and tolerance ladders 1e-3 to 1e-9, with every certification decision unchanged and returned probabilities identical in end-to-end testing. Scored against 40-digit mpmath folds of the engine’s own operators, the compiled fold errs at worst 2.4e-14 where the BLAS chain it replaces errs at 8.8e-14: the shift is dominated by the rounding of the code being removed, and the new fold is three to seven times closer to exact. A numba-less install is exactly 0.0 from the old code.

[1.0.10] - 2026-09-05

Changed

• _expm_stack asks whether a stack is anti-Hermitian in one compiled pass rather than five full-stack temporaries and about seven traversals, and builds K = 1j*Om only after the branch is taken. The framing stage falls 4.3-5.3x; a whole call falls about 1.1x, the difference being that the framing is only 6-17% of a call. Probabilities are unchanged: exactly 0.0 difference across 128 configurations spanning both profiles, two to five flavours, four slab counts and all four Magnus orders.

  The kernel compares squared magnitudes and takes one sqrt at the end, so scale and dev can differ from NumPy’s by up to 2 ulp. Computing abs() per element instead is exact in intent but measured 0.47-0.78x – slower than the NumPy expression it would replace – and is not value-identical either, so there is no exact compiled alternative to prefer. Those two scalars are used only to choose between the identity branch, the anti-Hermitian branch and the scipy fallback; the branch could change only for an input whose deviation-to- scale ratio sits within 2 ulp of 1e-12, four decades from where real input lands. Neither value reaches a warning or a returned quantity.

  NaN and inf route exactly as before. A plain maximum loop would not have done that – it skips NaN where numpy.max propagates it, which would have sent a NaN stack to the eigensolver instead of the scipy fallback – so the kernel carries a running sum of squared magnitudes, which vectorizes and cannot cancel, and checks it once at the end.

[1.0.9] - 2026-09-05

Changed

• Four- and five-flavour exponentials go to a batched Jacobi eigensolver rather than numpy.linalg.eigh. supports_dim now answers yes for 2 through 5. The reasoning it used to carry – that a 4x4 or 5x5 Hermitian eigenproblem has no practical closed form, so those dimensions keep eigh – had a true premise and a conclusion that did not follow: what made them slow was never the missing closed form but eigh’s fixed per-matrix LAPACK overhead, about 2.3 us on a 4x4, two thirds of a whole d = 4 pass. The kernel is a cyclic complex-Hermitian Jacobi sweep that warm-starts each matrix from its predecessor’s eigenvectors – consecutive matrices are consecutive slabs of one energy, so they arrive nearly diagonal – and re-orthonormalizes that basis at every step, which is load-bearing rather than tidy: without it non-unitarity compounds along the chain, 2.3e-11 after 13 000 matrices against 3.9e-14 with it.

  Marginal cost per slab falls 1.81-1.95x at four flavours and 1.46-1.59x at five, end to end on both profiles, arms interleaved in one process.

  Unlike every other kernel in this series, this one is not bit-identical, and it is not meant to be: it replaces an iterative eigensolver with a different iterative eigensolver. Worst full-matrix probability shift measured 1.557e-12. It is held instead to eigh’s own accuracy class, and to the same 10x bar the suite already applies to the closed-form kernels: across 364 cells spanning the norm sweep from 1e-150 to 1e4, crossed with clustered spectra down to exactly degenerate, its error stays within 6.4x of eigh’s, both scored against scipy.linalg.expm. Two and three flavours are untouched – max|dP| is exactly 0.0 there, and their timings do not move.

  Being backward stable with no conditioning cliff, it needs no SEV_TOL gate of the kind the closed forms require; it returns sev = 0.0, escalating to eigh only if a matrix fails to converge in 30 sweeps, which no census has observed. It terminates when a sweep performs no rotation, at a per-element threshold of 1e-30 relative. That threshold is not the prototype’s: the original squared off-norm floor sat inside the rounding equilibrium, so long chains declined to eigh almost always and the speed-up vanished silently.

[1.0.8] - 2026-09-05

Changed

• The order-4 Gauss-Legendre Omega is built by a compiled kernel that makes one pass over each slab’s two node samples, fusing three steps the NumPy form pays separately: the constant-sample equality test, the commutator, and the linear combination. The two samples are strided views, so the NumPy route copies both before the commutator can read them; the kernel reads them where they lie, and the equality scan stops at the first differing element rather than comparing every one. Marginal cost per slab falls 1.52x at two flavours and 1.23-1.36x at three, less at higher flavour counts where the eigendecomposition dominates. Orders 2, 6 and 8 are untouched and keep the commutator kernel of 1.0.7.

  Unlike 1.0.6 and 1.0.7, this one introduces no new divergence between a numba-less install and a numba one: the kernel reproduces the compiled commutator’s accumulation order and the same association of the scalar factors, so its output is bit-identical to the expression it replaces – exactly 0.0 difference across every shape, flavour, width form and edge case tested. Anything the kernel cannot take – another dtype, another node count – falls through to that expression unchanged.

• Orders 6 and 8 get the same treatment. Order 6 collapses about nineteen full-stack temporaries into eight d x d scratch buffers hoisted out of the slab loop; order 8, about sixty-five into nineteen. NumPy evaluates each operation across every slab before starting the next, so each intermediate is a full array written to memory and read back, while the kernel finishes one slab before moving on. Marginal cost per slab falls 1.29x at order 6 and 1.75x at order 8 end to end, median over eight profile/flavour cells. Both are bit-identical to the expressions they replace, and the six chained commutators of order 8 are each stored where NumPy stores one rather than algebraically re-folded.

  Two rounding differences had to be worked around to reach that, both worth knowing if these kernels are ever extended. NumPy divides a complex array by a real scalar through Smith’s algorithm, which multiplies by the reciprocal; numba divides componentwise, and about two thirds of two million random values differ in the last bit. Every X/c in these kernels is therefore written X*(1.0/c). Separately, numba’s ** does not round as numpy.float64.__pow__ does – 26% of random doubles differ in the last bit at **3 – so the order-8 weight-times-power scalars are precomputed as module constants in Python rather than raised inside the kernel.

  Order 2 was prototyped, measured and deliberately left alone: its branch is a single ufunc call with no boundaries to fuse, and at two flavours a kernel loses outright (0.82-1.04x per call), the dispatch’s own fixed cost exceeding what streaming saves on a branch costing 0.01-0.05 us/slab.

• The set of accepted pass-through keyword names is computed once and cached rather than rebuilt by inspect.signature on every public call. It was being rebuilt twice per call, once by the entry point and once by the osc_prob it delegates to, at about 221 us each; a single-point call is that much cheaper. This is fixed overhead, so it is invisible on large grids and worth most to callers making many small calls. cache_clear remains available for tests that add or remove parameters at runtime.

[1.0.7] - 2026-09-05

Changed

• The Magnus term recursion and the Gauss-Legendre schemes compute their commutators in a compiled kernel that fuses both matrix products into one accumulation. magnus.commutator itself is unchanged and stays the general pure-NumPy form; the kernel is a private path, and anything it cannot take falls through to the same expression. Marginal cost per slab falls 3.1-3.4x at two flavours and 2.0-2.2x at three at order 4, less at higher flavour counts, and not at all on the cumulative path.

  As with 1.0.6, a numba-less install is no longer bit-identical to a numba one on these paths: worst observed shift 6.7e-14 across 36 configurations, every refinement decision unchanged.

[1.0.6] - 2026-09-05

Added

• Order-8 Gauss-Legendre collocation integrator, on four nodes. 'gl' now reaches orders 2, 4, 6 and 8 from 1, 2, 3 and 4 Hamiltonian evaluations per slab; MAGNUS_EXP_ORDER_MAX_GL rises from 6 to 8, and a request above 8 raises as before. Orders 6 and 8 use the commutator-optimal forms of Blanes, Casas & Ros, BIT 42, 262 (2002), needing three and six commutators – the fewest possible at each order. Verified against a DOP853 reference: local error converges as h^9 (measured slopes 8.3-9.7 over three Hamiltonians), and on a PREM chord the row sums hold to 1.8e-15.

Changed

• The Gauss-Legendre schemes are described as collocation integrators throughout. Several places had called them “commutator-free”, which is the distinct family of Blanes & Moan (2006) and Alvermann & Fehske (2011) that replaces commutators with products of exponentials; these schemes are built from commutators, three at order 6 and six at order 8.

• The separable energy scan composes its slab operators in a compiled kernel rather than a Python loop. The association is unchanged – the same left fold, earliest slab rightmost – so this is not a reassociation; the kernel accumulates each matrix element as a compiled scalar sum where BLAS orders the same arithmetic its own way. Marginal cost per slab falls about 3.5x at three flavours and 2.3x at two; the gain at four and five flavours is smaller and was not resolved on the machine used. Without numba, or on a dtype the kernel was not built for, the original loop runs.

  One consequence is worth knowing: a numba-less install was bit-identical to a numba install on this path and may now differ at the 1e-14 level – worst observed 1.28e-14 across 16 scan configurations, with every refinement decision, warning and slab count unchanged. See docs/source/performance.rst.

[1.0.0] - 2026-08-11

Fixed

• The sterile states felt no matter at all on the NSI route. osc_prob_matter_nsi built the standard part of its matter Hamiltonian as a literal diag([1, 0, 0, 0]), so the sterile states got zero where they carry -V_NC = (r/2) V_CC. With every NSI coupling set to zero the route has to reproduce osc_prob_matter_std_potential exactly, and did not: it differed by 5.2e-02 at four flavors and 5.1e-02 at five. Three flavors agreed all along, which is why nothing caught it – the only test that can see it needs a fourth state. Anyone who has published 3+1 or 3+2 NSI-in-matter numbers from an earlier build should re-run them.

• The solar LIV routes read the Sun’s electron density as a mass density. osc_prob_Nnu_sun_liv builds its profile from NUM_DENSITY_E_SUN_CENTRAL, an electron number density, but forwarded a density_is_of_number_of_electrons flag defaulting to False. Against osc_prob_Nnu_sun with the couplings zeroed they differed by 0.69, 0.69, 0.45 and 0.43 at 2, 3, 4 and 5 flavors. osc_prob_Nnu_sun and ..._sun_nsi never exposed those flags, which is why only the LIV family was affected.

• A negative energy returned the antineutrino probability. E < 0 flips the sign of the whole Hamiltonian, which is CP conjugation, so the call returned a unitary, entirely plausible answer to a question the caller did not ask – matching nubar=True to 1e-15 and differing from the intended answer by 2.3e-02 in vacuum and 4.4e-02 in matter. Non-positive energies are now rejected.

• A NaN density was reported as a units mistake, because nan == 0.0 and nan >= threshold are both False and the guard fell through to its warning. Non-finite densities are now rejected, naming the real problem.

• Notebook 12’s four- and five-flavor ground truth built its comparison Hamiltonian with the same zero-sterile literal, so solve_ivp was integrating the wrong problem and the error column blamed the strategy for the reference. The four-flavor standard case goes from err_magnus = 1.73e-04 to 2.19e-06.

• The quickstart’s and README’s code examples could not be run – one used names the page defines fifteen lines later, others passed ... as a Hamiltonian body. Both are now executed by the suite.

• The documented solar averaging reduction outlived the script it came from. The averaged-probability page, notebooks 13, 14 and 23 and the tutorials index all asserted “instantaneous 1.4e-03 -> averaged 2.6e-05, a 53x reduction”, citing adversarial_batteries/avg_check.py. That script had gone stale at 0bf3a40, which is on main and predates this audit; re-run, it gives 6.000e-04 -> 7.110e-04 on the log-linear BS05 ray notebook 13 actually uses (8.889e-04 -> 6.051e-04 on the cubic-spline variant). All four rows of the table are now measured values.

  The rule stated alongside them was the worse problem and is withdrawn: the docs told the reader to read the ratio of the two columns as a diagnostic, errors shrinking “more than twentyfold” being phase and the rest envelope. A finite-window mean is an estimator with a bias of its own, and on a profile whose density varies across the window that bias does not shrink as the window widens – on the solar ray the window mean drifts 0.5924 -> 0.6023 between six and forty-eight oscillation lengths, away from a limit rather than toward one. So the ratio is meaningful only at fixed matter conditions, as in notebook 23, and the page now says so. To obtain the averaged probability, ask for it: average=True evaluates the decohered limit in closed form and reproduces the adiabatic MSW expression to 3.33e-16 across 1–20 MeV.

• The front page’s solar timing comparison outlived its measurement too. The README and the docs index both read “40 averaged energies in 0.66 s, against 131 s for 12 instantaneous ones from nuSQuIDS”. That was measured, but on a superseded design of notebook 25’s section 10; the frozen external_solar_nusquids.json now sweeps four solver tolerances over a 200-point grid targeting 40 energies, and contains no 12-energy run and nothing at 131 s. Both pages now say what the notebook says and what its data supports: nuSQuIDS needs about ten minutes merely to reach the tolerance at which its output is a probability at all (568 s, the cheapest setting the generator marks physical), and a further factor of N to average the phase away. Magnus’s own 0.66 s is a live timing and re-measures at 0.68 s, which is run-to-run noise on one machine and is left as it stands.

• osc_prob_vacuum documented twelve refinement knobs as forwarded that it deliberately drops. t_slab_edges, magnus_exp_order, integration_method, rtol, atol, both growth_factor_*, max_num_loops and the four min/max_n_* bounds each said “Forwarded to osc_prob_energy_baseline /osc_prob”, while the body forwards only n_jobs, validate_input, verbose, save_log and file_log. Not forwarding them is correct and the code says why – a vacuum Hamiltonian is constant in position, so every point is exact with a single slab and there is nothing to refine – so the docstring was the wrong text, and it is the one a caller reads before passing rtol=1e-9 and getting silence. No behavior changes. filename_log and close_file_log_upon_exit were documented the same way and are likewise never read.

  Found by generalizing tests/test_angles.py’s “declares the parameter, never reads it” AST check from angles to every parameter in src/magnus. The other 51 hits are the legitimate interface-signature pattern – _td builders that take l and ignore it, constant-density closures – or already documented as inert, such as hamiltonian_2nu_liv’s nubar, which has nothing to conjugate because the 2-flavor LIV rotation carries no phase.

• SEV_TOL’s documented accuracy guarantee was not the measured one. The docstring claimed 2e-13 absolute across everything the gate admits; re-measured on two spectrum families with 40 random bases per rung, the admitted range reaches 2.0e-12, and even the m <= 1.1e3 corner the claim was calibrated on exceeds 2e-13 in about 1% of bases. No value of the constant could have rescued it, because the test that pins the window requires the gate to sit above the very cell where the claim already fails, so the number was corrected instead: 5e-12 across the admitted range and 5e-13 in that corner, each about twice the worst measured. The gate itself is unchanged at 1e4 and no result moves.

  Two earlier calibrations of this constant appeared to contradict each other – one putting the first unsafe cell at m = 1.1e5, the other at 4.4e3. They do not: they used different spectrum families under the same “scale 1e2” label, [-s, -s(1-d), s] spanning 2s against [0, d, S] spanning S, a factor of four in m. Compared at equal m they agree. m = tr(X^2)/6 is a spectral invariant and is the only fair label; “scale” is not.

• The sterile matter term and the Earth’s density describe different media. L37 made the Earth’s Y_e a function of radius and derived r = (1 - Y_e)/Y_e from it layer by layer, because the two are one statement about composition. The sterile states’ entry in the matter projector, r/2, could not follow: it is a single matrix for the whole chord, so it kept taking the caller’s scalar, which defaults to 1.0 – isoscalar matter, i.e. Y_e = 0.5, precisely the uniform composition L37 replaced.

  So on every Earth chord with a sterile state the density and the projector disagree by construction. Measured at costhz = -0.95 with s14 = 0.15, s24 = 0.10 and D41 = 1 eV^2, the isoscalar projector differs from one built with the core’s own r = 1.1478 by 2.1e-02 in P(numu -> numu) – twenty times the default tolerance, and silent. Three flavors never saw it, because the projector’s sterile block is empty; that is the same reason A2b survived a max-effort review.

  No single r is right for a chord that crosses iron and rock, so this is reported rather than resolved: the Earth wrappers now raise globaldefs.SterileMatterCompositionWarning when the two disagree by more than 2%, and the message names the path-averaged ratio for the chord that was asked for, which silences it. electron_fraction=0.5 with the default 1.0 makes them agree exactly and reproduces the older uniform composition. No number changes unless you act on the warning.

  The default was deliberately left at 1.0 rather than moved to the path-averaged value, because no single r is right for a chord crossing iron and rock and a self-consistent looking default would hide that. The real fix is a position-dependent projector, H_matt(l) = V_CC(l) x P(l), which is a structural change to the scalar-times-constant- matrix factorisation the matter path is built on; it is tracked for a future release as issue #47.

  The docstrings were also wrong about this. Twelve Earth wrappers said the ratio “must match the value given to vcc_func_from_rho_func”, which on the Earth path is unsatisfiable – that function is handed a per-layer ratio the caller never supplies.

• The sterile states’ matter entry was unreachable on the Sun, and locked isoscalar. osc_prob_{4,5}nu_sun and their _nsi/_liv variants never exposed ratio_number_neutrons_to_protons, and delegated without forwarding it – two of them passing a hardcoded 1.0. So r/2, the sterile entry in the matter projector, was fixed at isoscalar for a medium that is nothing of the kind: the Sun is hydrogen-rich, Y_e = (1 + X)/2 runs from about 0.68 at the center to 0.88 near the surface, and r = (1 - Y_e)/Y_e from about 0.47 down to 0.14. Unlike the Earth, where 1.0 at least sits among the layer values, for the Sun it is outside the physical range entirely.

  The six wrappers now take it and forward it. Left at the default the averaged solar survival probability moves by about 4e-03 at s14 = 0.4, above the 1e-3 default tolerance; at the nominal s14 = 0.15 it is 7.9e-04, just under. The solar profile is a fit to the electron number density, so Y_e is already inside it and there is nothing for the library to derive r from – it has to be stated, which is why this is exposed rather than defaulted to a solar value. No number changes unless you set it.

• t_breakpoints was unusable on every Earth wrapper. Those wrappers place slab edges on the PREM shell crossings themselves, so an argument of the same name arrived in **kwargs and collided: the caller got got multiple values for keyword argument 't_breakpoints' raised two layers down. The keyword is listed as forwardable in this package’s own unrecognized-keyword message, so it was reachable and broken. A caller’s breakpoints are now merged with the PREM crossings rather than either replacing the other – dropping the crossings silently would be the very defect t_breakpoints exists to prevent. Pass t_slab_edges to place every edge yourself.

• n_jobs is documented as not being a pure performance knob. Splitting slabs across workers changes the order the arithmetic is done in, and the refinement ladder’s stopping test compares successive levels, so it can stop one level earlier or later than the serial run. Measured on a 3nu PREM chord over eight energies, serial against two workers differs by 1.2e-03 at the default rtol = 1e-3, 6.6e-08 at rtol = 1e-6 and 5.6e-11 at rtol = 1e-9 – within the tolerance asked for, but not bitwise, which the docstring previously gave no hint of.

• The command line could not reach half its own parameter sets. --osc-params-set listed its choices by hand and offered only the NuFit 6.0 entries, so after the default moved to 6.1 anyone asking for inverted ordering on the command line silently dropped a release behind the default. The choices now come from globaldefs.OSC_PARAMS_PREDEFINED itself.

Added

• An angles keyword on every function that takes a mixing angle. It selects the convention the angles are stated in: 'sin' (the default, and the only behavior before this) their sines, 'sin2' their sines squared – the form global fits are published in – 'rad' the angles in radians, or 'deg' in degrees. Under 'deg' the CP phases are read as degrees too; under the other three they stay in radians, a sine being no way to state a phase. Ninety-five functions take it, across hamiltonians, oscprob, oscprobstd, the magnus prob command line and globaldefs.load_nufit_params, and the four routes agree to 2.2e-13 end to end. The default is a pass-through, so nothing that does not ask for it changes.

  It exists so a published parameter set can be typed in as published: s12=0.308, ..., angles='sin2', or s12=33.76, ..., dCP=212.0, angles='deg', rather than square-rooted and degree-converted by hand at the call site. It also matches the keyword NuOscProbExact uses, so the two codes can be driven from one parameter set when they are compared.

  Four guards come with it, because a mis-stated convention is otherwise silent: an unrecognized value, a sine outside [-1, 1], a negative 'sin2', and an angle above 2*pi under 'rad' all raise; a whole parameter set under one degree with 'deg' raises globaldefs.MixingAngleConventionWarning, since theta_13 at about 8.5 degrees is the smallest angle anyone measures and values that small are sines. 'rad' and 'sin' are the one pairing no bound can separate – theta_12 = 0.589 against sin theta_12 = 0.556 – which is why load_nufit_params takes the keyword too: convert once, at the source, rather than by hand between the two.

• matter.matter_potential_projector is exported. It is the one definition of the matter term’s structure, and every place that rebuilt that structure by hand instead got the sterile entry wrong. While it was unexported, autoapi did not document it and references to it from public docstrings resolved to nothing.

• Input validation for the engine’s own knobs: non-positive n_slabs, min_n_slabs, n_tpts_per_slab, and a floor set above its own ceiling, were all accepted and then quietly ignored, so a typo looked like a setting that had been honored. rtol and atol were already guarded; these now match.

• A PEP 561 py.typed marker, so the annotations the public API already carries are visible to mypy and pyright instead of being discarded.

• CITATION.cff, so GitHub offers “Cite this repository” and reference managers can read the metadata.

Changed

• The Earth’s electron fraction is now resolved per PREM layer, and results change. Every Earth entry point assumed Y_e = 0.5 – exactly isoscalar matter, which nothing in the Earth is. PREM is a density model and carries no composition, so Y_e = <Z/A> has to be supplied; it is now taken per layer, at radii that are already PREM boundaries:

  layer	radii [km]	Y_e	material
  core	r <= 3480	0.4656	iron
  mantle	3480 - 6346.6	0.4957	peridotite
  crust	6346.6 - 6368	0.4952	granitic
  ocean	r > 6368	0.5551	seawater (H has Z/A = 1)

  The correction tracks how much core a chord crosses: at 1 GeV, P(nu_mu -> nu_e) falls to 23% of its previous value at cos(theta_z) = -1, and moves by about 1% at -0.4. Anyone with published Earth numbers from an earlier build should re-run them; passing electron_fraction=0.5 reproduces the old uniform composition exactly.

  Each layer is settable (electron_fraction_core and friends). Combining a per-layer value with the uniform electron_fraction is refused rather than silently resolved, and Y_e is validated as a fraction in (0, 1] – 0.0 and 5.0 used to be accepted, returning answers 0.51 and 0.74 away from the default. The neutron-to-proton ratio is derived from Y_e in the density conversion, r = (1 - Y_e)/Y_e, so a caller can no longer describe an iron core with isoscalar neutrons.

  Two caveats are documented rather than guessed at: the crust value differs from the mantle by 0.1%, so it exists for explicitness rather than effect; and PREM’s ocean is a global average that a land-based baseline does not cross, for which electron_fraction_ocean=Y_E_CRUST_PREM is the right setting.

• The solar LIV routes no longer accept four parameters they ignored. osc_prob_Nnu_sun_liv took electron_fraction, ratio_number_neutrons_to_protons and two density flags and used none of them: passing electron_fraction=0.25 changed the answer by exactly zero. The solar profile is the standard exponential fit to the electron number density, so the mass-density conversion those describe never runs and Y_e is already inside the fit. osc_prob_Nnu_sun and ..._sun_nsi never exposed them.

• The default oscillation parameters are NuFit 6.1, and there is now only one set of them. Omitting oscillation parameters used to fall back to a second copy of the numbers built from NuFit 6.0 constants, while load_nufit_params() with no arguments returned 6.1 – so the same script got different answers depending on which door it came through, by 4.0e-03 in probability at 1 GeV over 1300 km. OSC_PARAMS_DEFAULT is now derived from load_nufit_params, so the two cannot drift apart again. OSC_PARAMS_NU_FIT_6_0_SK_NO and ..._SK_IO remain available by name for anyone reproducing an earlier number. Results that relied on the implicit default will change.

• The publish workflow gates on the tests, on twine check --strict, and on the release tag matching the packaged version. It previously built and uploaded on a release with no check of any kind, and PyPI does not allow re-uploading a version.

• CI tests Python 3.13 as well as 3.10-3.12, since requires-python has no upper bound and pip will install on it either way.

• Notebook 27, nine animated scenes, with the clips committed. Four are the scenes NuOscProbExact’s notebook 19 draws, computed here so the two can be read side by side; the other five need something a closed-form slab code does not have – a refinement ladder deciding it has converged, a front that travels, an observable that is an average rather than a value, and a Hamiltonian that varies along the path. The notebook draws stills by default and hides the rendering behind RENDER = True, so CI never pays the hour it costs. Seven shrunk GIFs are tracked in img/ (14.7 MB); the raw renders go to img/raw/, which is gitignored, so a later render cannot silently replace 14.7 MB of committed files with the 224 MB it produces. tools/make_demo_video.py owns the encoding.

• Notebook 25 became an arbiter of when each code wins, rather than a list of timings: reach (a slab product’s error floors at 2.5e-11 on a smooth profile and then rises, while the Magnus expansion continues to 2.9e-13), generality (five flavors, where there is no comparison to draw), and pre-packaged observables (a solar average in 0.66 s against 130 s). With it, supernova-shock comparisons against other codes, 3+1 and NSI in both the solar and shock settings, and probability-vs-energy panels for the shock.

• A README image gallery, ten figures lifted from the executed notebooks by extract_gallery(), so the front page shows the answers rather than describing them and cannot drift from what the notebooks produce.

• An expansion-order section in notebook 24: what the truncation order buys in accuracy, and what it costs in correctness, which is nothing – every truncation lives in the Lie algebra, so the operator is unitary exactly rather than to the accuracy of the truncation. With it, which engine the dispatcher picks and why, measured across four decades of tolerance.

• t_breakpoints, n_slabs and cumulative reach the BSM wrappers. The keywords that decide a hard profile were in no signature that a caller of osc_prob_matter_nsi or osc_prob_liv could see, so the comparison those wrappers exist for could not be made on equal terms.

• A 'constant' engine: a position-independent Hamiltonian is answered in one batched exponential instead of one osc_prob call per point. When the matter potential does not vary with position, the Magnus series terminates at its first term — Ω₁ = −iHΔ and every higher Ω is a nested commutator of H with itself, hence zero — so U = exp(-iHΔ) is the exact answer and a whole energy scan is one stacked exponential.

  This case was previously turned away on purpose: _osc_prob_scan_separable_dispatch bailed on not isinstance(VCC_func, Callable), its docstring saying “a constant potential falls back to the generic path”. So the easiest Hamiltonian there is took the slowest route available — a 60-energy scan made 18,000 osc_prob calls per 300 repetitions, each rediscovering the same constancy and paying the full wrapper and refinement-ladder overhead.

  Measured against the route it replaces, interleaved with a control that came back at 1.00×:

  flavors	matter scan	vacuum scan	single point
  2ν	17.3×	24.7×	2.0×
  3ν	15.5×	18.9×	2.1×
  4ν	7.2×	7.4×	1.4×
  5ν	6.0×	6.2×	1.4×

  4ν and 5ν gain less because they are on eigh rather than the Cayley–Hamilton kernel, which covers dimensions 2 and 3 only. A 3ν constant-density scan is now 1.10 µs per energy against NuOscProbExact’s 1.44 µs batched and 13.25 µs looped; a single point is 33.8 µs against its 19.9 µs, the remainder being wrapper parameter resolution rather than arithmetic. Results are bit-identical to the per-point route on every flavor count and both neutrino signs, and n_slabs, n_tpts_per_slab, t_breakpoints and rtol/atol are accepted and ignored because they can only ask for refinement of something already exact.

  PREM and exponential profiles are untouched and keep separable/magnus /hybrid: their potential varies with position, and a constant-H engine that captured one would propagate the whole trajectory with a single exponential — wrong by O(1) while still perfectly unitary. A test asserts the engine identity, not merely the numbers.

• A compiled Cayley–Hamilton backend for the matrix exponential, selected by magnus.magnus.EXPM_BACKEND. np.linalg.eigh costs ~1.27 µs per 3×3 whatever the stack size (measured 1.268 µs at N=108, 1.279 µs at N=4096 — flat, because it loops over LAPACK internally instead of vectorizing). The new magnus.expmkernels applies to K the polynomial interpolating exp(-iλ) on its spectrum instead: no eigenvectors, and the eigenvalues in closed form. 6.8× on the exponential at N=108 (162.6 → 23.8 µs), 7.3× at d=2, and 2.11× end to end on a 60-energy PREM scan (9291 → 4409 µs, i.e. 73.5 µs per energy).

  The gap between 6.8× and 2.11× is Amdahl’s law: the exponential is about a third of a slab pass. Quoting the 6.8× as a package speed-up would be quoting the wrong number.

  'auto' (the default) uses the kernel for 2×2 and 3×3 when numba is installed and eigh otherwise, and cannot fail; 'numba' makes a missing numba an error rather than a silent downgrade; 'eigh' is the reference route. Dimensions 4 and 5 keep eigh — there is no practical closed form for a 4×4 or 5×5 Hermitian eigenproblem, so 4ν and 5ν stay correct and are not accelerated. numba is an optional dependency (pip install 'magnuspy[fast]'), costing ~90 ms of import magnus when present plus a one-off ~0.7 s compile per kernel, cached to disk thereafter.

  Switching backend moves probabilities by at most 4.6e-15 across PREM chords, energy scans, NSI resonances, constant density and vacuum.

  Degeneracy is the whole risk in such a scheme, and two facts remove it. A Hermitian matrix is never defective, so matching exp on the distinct eigenvalues is already exact and the confluent (Hermite) form is not needed. And with eigenvalues sorted and the spectrum shifted to put the median at zero, the one ill-conditioned coefficient multiplies a matrix whose norm shrinks with the same gap, so its contribution is bounded by ε·gap and vanishes as the gap closes. There is therefore no tolerance, no crossover, and no near-degenerate branch to place: the error is 1e-16 at splittings of 1e-2, 1e-6, 1e-10, 1e-14 and exactly zero alike.

  Two things this cost, both now pinned by tests. The closed-form eigenvalues degrade to ~4e-9 at an exact degeneracy (arccos has infinite derivative where a repeated root sits) and the exponential stays at 2.5e-16 anyway, because interpolation error is second order in the displacement of a coalescing node; both halves are asserted, the sloppy one included. And the kernel must read the lower triangle, because that is the one eigh reads (UPLO defaults to 'L') and _expm_stack admits input anti-Hermitian only to 1e-12 — a kernel reading the upper triangle exponentiates a different matrix on such input and the two backends diverge by ~2e-12, large enough to matter and small enough to look like rounding.

• Palindromic density profiles are exploited on Earth chords. A chord through a spherically symmetric Earth meets every radius twice, so its density profile reads the same from either end. The Magnus core now evaluates the Hamiltonian on the first half of such a slab chain and derives the rest by reversal, halving the calls to the caller’s H_func. Worth 1.4x-1.67x on a single point and 1.56x-1.64x on an energy scan when that Hamiltonian is expensive; plain PREM, whose density lookup is cheap, pays about 10% for it. New public magnus.magnus.USE_PALINDROME (module switch, True) turns it off, and magnus.magnus.palindromic() is the predicate.

  The saving is halved Hamiltonian evaluations and nothing else, so it is worth what that Hamiltonian costs. Standard PREM scans are unaffected: they are answered by the separable engine, which already evaluates the profile once and shares it across energies – the same saving, taken earlier.

  Symmetry is declared by the Earth entry points, where it is a fact of chord geometry, not detected: detecting it would need the very evaluations the optimization skips. There is deliberately no user-facing way to declare it of an arbitrary profile.

  This moves Earth single-point results by up to 8.6e-15 relative. The mirrored slab’s nodes are reached by a different floating-point expression for the same real number, so the change is inherent rather than incidental. USE_PALINDROME = False reproduces the previous numbers exactly.

Changed

• Per-call overhead cut across the wrappers, by caching what is pure and cheapening what is common. The largest single item in a single-point profile was hamiltonian_3nu_vacuum_energy_independent at ~15 µs, rebuilding the same PMNS matrix on every call; it is a pure function of eight scalars and is now memoized, handing back a copy so a caller writing into the result cannot poison the cache. Likewise the constant-density branches of matter.vcc_func_from_rho_func (the callable branches are deliberately not cached: they return a closure over rho_func that callers tag with is_exp_density_profile, so caching those would trade microseconds for an aliasing bug). _n_required_params cached inspect.signature weakly against the function — 42 µs of cumulative time per call, and once per point on the routes that legitimately loop. isinstance(x, typing.Callable), which routes through ABCMeta.__instancecheck__, replaced by the callable() builtin at 23 sites (~9 typing.__subclasscheck__ calls per invocation), and scalar fast paths added to _normalize_energy_L and the density-units guard. No behavior changed by this entry.

• A verbose run (verbose >= 1) takes the per-point route. The banner and run-parameter dump describe quantities the batched engines do not have — magnus_exp_order, slab counts, tolerances — so emitting them from a batched path would report a refinement ladder that never ran.

• rtol/atol are documented for what they are: a stopping criterion, not an accuracy guarantee. The ladder halts when two successive levels agree; it never estimates the error of the answer it returns, which is a weaker promise than a stepping ODE integrator’s rtol makes. Corrected in osc_prob, adiabatic.hybrid_propagator, the CLI’s --rtol/ --atol help, README.md, architecture.rst (which said “until rtol/atol is met”), and a new section of implementation_details.rst that the others link to. No behavior changed by this entry.

• convergence_info reports what the ladder did. Alongside the existing n_slabs/n_tpts_per_slab it now carries n_slab_edges and n_slab_edges_previous (which make the real refinement step visible), n_slabs_previous, n_tpts_per_slab_previous, last_gap (None when only one level was ever computed), n_agreements, and tolerance_achieved – the programmatic form of ToleranceNotAchievedWarning.

  It deliberately carries no error estimate. Converting the gap into one by Richardson extrapolation, as the sibling NuOscProbExact does, was measured and rejected: Magnus has no stable convergence order (fitted on Earth chords it scatters from 1.4 to 7.2 against nominal orders of 2 and 4), and because breakpoints make the effective refinement ratio as low as 1.06 rather than 1.5, dividing by r^p - 1 under-reports the true error by 6-20x even where the power law holds exactly. Under-reporting is the dangerous direction.

• oscprob.BATCH_WORKING_ENTRIES lowered from 4,194,304 to 65,536 (about 67 MB to about 1 MB). The batched scan engines are memory-bound, and the previous value was large enough that their working set spilled cache. Measured across fifteen workloads on three engines, the new value is 1.19x-1.38x quicker on Earth energy scans, 1.06x-1.16x on cumulative baseline scans, flat on short scans and on the interaction-picture engine, and never slower anywhere. Bit-identical at every budget tested – tiles are independent and only concatenated – so this changes no result. Peak memory of a long scan drops accordingly.

Fixed

• The memory guard read the host’s free memory, not the cgroup’s. Inside a container the two are unrelated, so a request that would be killed by the cgroup limit was waved through by a guard reading a number that did not apply to it. _available_memory_bytes now takes the minimum of the two, walking the cgroup ancestor chain and handling the v1 “no limit” sentinel.

• The sterile state felt no medium, in six places. Four inline copies of the matter projector, plus the two flavor-specific builders, wrote the 3+1 matter term by hand and gave the sterile state a zero where it carries -V_NC. All six now come from matter.matter_potential_projector, and a test fails if a seventh copy appears. On a PREM chord this was worth 0.29 in probability, and it was flat in the requested tolerance, so no amount of refinement revealed it.

• The evaluation-mode cache answered a question it was never asked, and the guarded version of it was dead code. probe_eval_mode’s 'constant' verdict means “sampling A across [t0, t1] gave the same matrix every time” — a property of the function on an interval, not of the function. Keyed on the function alone, a mode learned on a short baseline was served for a long one, and _evaluate_A’s 'constant' branch then broadcast a single sample over a profile that genuinely varies, with no spot-check (unlike the 'vector'/'scalar' hints, which self-heal).

  A two-layer profile written with the natural short-circuit — if every requested position falls in one layer, return that layer’s matrix — probes as 'constant' on a short interval and 'vector' on a long one. Looping baselines shortest-first then gave P_ee = 0.906249 against a correct 0.903424, row sums exactly 1.0, no warning, and the right answer if the loop ran the other way round. Now keyed on (function, t0, t1): measured loop-order dependence 5.8e-02 → 0.00e+00.

  Three further defects were the same defect. oscprob._eval_mode_for was a second, unguarded copy of this cache that reached across a module boundary into magnus._EVAL_MODE_CACHE; it guarded its store but not its lookup, so a callable Hamiltonian that cannot be weakly referenced (__slots__) or hashed (a @dataclass, or anything defining __eq__) raised TypeError from inside osc_prob where it had worked before. The public magnus.cached_eval_mode, which was guarded, had zero live callers: its only call site sat in the not callable(H_func) arm of a ternary, and H_func is unconditionally rebound to a closure above it. The copy is deleted, the original does the work via a new key parameter, and the comment claiming the cache keys by identity is corrected — WeakKeyDictionary keys by the referent’s equality, which matters for a Hamiltonian defining __eq__.

  The interval key costs nothing measurable: a refinement ladder calling repeatedly at one interval still probes once, which is what the cache is for.

• A second max-effort review, run from a fresh session, found four more; all are fixed here. The first review below was written and verified by the agent that wrote the code, which is the reason to look again. None of the four changes a computed probability: they are a lost warning, an inconsistent error contract, an ignored backend request, and an unbounded cache.

  The DensityUnitWarning repair in the review below was half a fix, and its other half was a new defect. vcc_func_from_rho_func’s constant-density cache skips the conversion that the two unit guards live inside, and the earlier repair re-emitted only the density_matter_is_in_g_per_cm3=True arm. The arm it dropped is the dangerous one: an undeclared g cm^-3 density returns exactly the vacuum probability, so the warning is the only thing separating it from an answer. Re-emitting from the cache site cannot work either — warnings.warn’s stacklevel attributes the call to a different frame, and the frame is part of the interpreter’s registry key, so the imitation printed a second warning under the default filter where an uncached call printed one. A density that would trip either guard is now not cached at all, so both keep firing from where they always did. Measured against main: 3 of 3 identical calls warn under simplefilter('always'), 1 of 3 under the default filter, matching in both directions.

  The constant engine accepted three refinement parameters the ladder rejects. max_n_slabs=0, rtol/atol ≤ 0 and max_num_loops=0 were all answered here while osc_prob raised ValueError for each, so whether a bad parameter was reported depended on whether the density happened to be constant. The answers were never wrong — one exponential per point is exact either way — but the error contract was. _refinement_params_rejected mirrors osc_prob’s validation and declines, so the caller sees osc_prob’s own message rather than a second wording of it.

  EXPM_BACKEND did not cross a process boundary. loky re-imports magnus in every worker, where the switch is back at 'auto', so any call with n_jobs != 1 silently ignored it — and it is the only backend control the oscprob wrappers expose, the expm_backend parameter reaching no further than the Magnus layer. Worst for the one use the switch is documented for: a backend comparison run in parallel compared 'auto' against itself. Now carried by value into the worker and re-applied there.

  The evaluation-mode cache’s per-interval dictionary was unbounded. Weak keys bound the outer map, not the inner one, and a Hamiltonian defined at module scope never dies: 1000 distinct baselines through osc_prob retained 1000 entries, about 184 KB, for the life of the process. Now bounded at 256 and cleared wholesale, matching _VACUUM_H_CACHE and _VCC_CONST_CACHE — the ladder holds one entry no eviction can reach, and a scan uses each entry once, so neither population rewards a smarter policy.

  A refinement bound named the wrong parameter, and had since long before this branch. Found while mirroring osc_prob’s validation for the fix above: one of its two ceiling checks tested max_n_slabs while its message named max_n_tpts_per_slab. So max_n_tpts_per_slab was never validated at all — 0 and 1 were accepted — max_n_slabs was bounded at > 2 while the message three lines above promised > 1, and a caller who passed max_n_slabs=2 was refused in the name of a parameter they had not touched. The rule both messages encode is that each ceiling clears its own floor (min_n_slabs defaults to 1, min_n_tpts_per_slab to 2), so the condition was the wrong half: it now tests max_n_tpts_per_slab. max_n_slabs=2 starts working and max_n_tpts_per_slab <= 2 starts raising; nothing in the tests, the notebooks or the docs passed a value that newly raises.

  The docs gate had been red since 2debd51, and the branch did not know it. A :func: role in implementation_details.rst pointed at magnus.magnus._expm_stack, a private name autoapi does not document, which fails -n -W. It went unnoticed because every check of that gate on this branch was an incremental build, and Sphinx does not re-read a file it believes unchanged — so the warning could not reappear once its file had been cached. CI builds from a clean tree and would have caught it on the first push. Now double backticks, which is what the same file does for _osc_prob_scan_constant_h and for _expm_stack itself nine lines later. The gate is re-verified from make clean: 0 warnings.

  Examined and left alone: the Cayley-Hamilton algebra (det X, the divided differences, the root ordering and Z² re-derived independently against the code); the d = 2 kernel, which has no SEV_TOL gate but beats eigh two orders of magnitude beyond its documented range (1.9e-09 against 3.3e-09 at ‖K‖ = 1e7); float(VCC_func) on an array-valued density, unreachable because validate_input_battery rejects it first; and verbose=None, which already raised on main. The constant engine was re-checked against scipy.linalg.expm rather than against another magnus engine — the two in-package routes share _expm_stack and agree to exactly 0.0, which is no evidence at all — and matches to 3.4e-15 across ν/ν̄, three densities, per-point baselines, L0 ≠ 0, a single scalar point and vacuum.

• A max-effort code review of this branch found fifteen defects; all are fixed here or in the two commits below. Eight independent finder angles, every finding confirmed by execution rather than reading. The two that mattered most were invisible to the tests that were supposed to catch them.

  The compiled exponential was up to 7440x less accurate than eigh where a clustered spectrum meets a large norm — 2.7e-07 against 3.0e-11 at ||K|| = 1e5, because arccos has infinite derivative at u = ±1. It had been verified against random spectra at many norms, and separately at many eigenvalue separations at norm ~1; the damage needs both at once, which no single-axis sweep visits, and this file previously claimed the kernel was “the same order or slightly better at every norm” on that evidence. expmkernels.SEV_TOL now hands such matrices to eigh: worst absolute error over the whole separation-by-scale grid 8.7e-14, matrices declined on real work 0.00%, speed unchanged. The grid is now a test.

  Two tests were vacuous. The antineutrino-sign test — written specifically to catch the h_matt bug in the commit below — passed identically with the engine under test disabled. The agreement test compared the new engine against osc_prob’s other constant shortcut rather than the refinement ladder, so its tolerance could never fire; instrumented, the slab machinery ran zero times on either side. Both now assert the route they claim to compare, and both were verified to fail when their target bug is reintroduced.

  Also: expm_herm_stack ignored supports_dim and handed 4x4 input to the 3x3 kernel (error 2.4, unitarity 11.3, uninitialized eigenvalues) and segfaulted at d=1; the constant engine answered L < L0 with the transpose of the right answer (29% off, row sums exactly 1) where main raised, accepted magnus_exp_order=0, and bypassed the output-size memory guard; verbose=1 lost all its output; array-valued and 0-d-array parameters began raising in two new caches; and a cache hit silenced DensityUnitWarning after the first call.

  Four further findings inherited from earlier commits on this branch are recorded in docs/dev/HANDOVER_OVERHEAD.md and deliberately left for separate work.

• h_matt meant two different things depending on the potential, and the new engine walked into it. osc_prob_matter_nsi and osc_prob_liv rebound h_matt to VCC_func*h_matt on their constant-potential branch, then passed that name to dispatchers documented to take h_matt as “the constant matrix multiplying VCC_func(l)” — which multiply by VCC themselves. The separable engine never noticed, because the rebinding only happens on the constant branch it used to decline outright.

  The failure was invisible in the two ways that matter: VCC² is ~1e-25 rather than ~1e-13, so the matter term all but vanished and the answer stayed a plausible, unitary, nearly-vacuum probability; and VCC² has no sign, so the neutrino and antineutrino results came back bit-identical. Standard constant-density matter was unaffected (that call site passes the bare projector), so testing the headline case alone would have missed it. The scaled matrix now has its own name, and tests/test_engines.py compares the two routes across every scenario wrapper and both signs — verified to fail when the bug is reintroduced.

• A new engine was invisible to the cross-check. _CROSS_CHECK_FORCING’s forbid lists and ENGINE_FAMILIES did not know about 'constant', and since it answers before osc_prob_energy_baseline — which is what records the payload the independent expm reference is built from — enabling it silently removed the only non-Magnus oracle in the table. It is now listed in every other row’s forbid set, and shares the 'exact' family with expm rather than standing alone: the two use different exponential implementations but share the assumption that H is position-independent, and that assumption is the thing that could be wrong.

• _expm_stack’s docstring claimed the eigh route was exactly unitary, and it is not. U†U - I measures 4e-16 for a single 3×3 and 4e-15 for a stack of 4096 — growing with stack size, never zero, because reconstruction from eigenvectors rounds like any other floating-point product. The claim that probabilities “sum to 1 by construction” was the part worth correcting: they sum to 1 to about 1e-15, which is worth relying on, by rounding rather than by construction. No behavior changed by this entry.

• The refinement ladder could stop while the answer was still outside the requested tolerance, and report success. t_breakpoints (the ~14 PREM layer crossings) are re-inserted into every level’s grid, so at small counts the nominal refinement and the real one are different things: a nominal 2 -> 3 slab step is a 16 -> 17 edge step, a 6% refinement rather than a 50% one. Two grids differing by 6% agree for reasons unrelated to having converged, and np.allclose read that as success.

  An agreement now only counts when the two levels compared were genuinely different grids – len(t_slab_edges) must grow by at least :data:oscprob.MIN_EFFECTIVE_REFINEMENT (1.25, measured). Over 120 Earth configurations (costhz -0.15 to -0.99, 0.5-8 GeV, three tolerances, scored against a reference verified converged to 1e-13): one silent violation, 2.1x outside the tolerance asked for, and five more that missed but warned – all six now zero.

  This changes Earth results, by making them more accurate: one of the twelve bitident rows moves, the Earth single point, by 2.3e-04. It is inert wherever there are no breakpoints (solar is bit-identical), and does not reach the separable engine that answers Earth scans, which was checked separately and does not have the defect. Cost is about 11% of wall clock (a 60-energy Earth scan goes from 9 ms to 10 ms) despite the median slab count rising from 9 to 21, because the added levels are the cheap small ones.

• The position-profile cache no longer hands out writable arrays. Values are returned by reference to every later caller asking for the same position grid, so a write through any one of them would have silently changed what the others received – and the cache sits under the matter term of the Hamiltonian, so the symptom would have been a wrong probability with nothing raised. Cached arrays are now marked read-only, turning that into an exception at the point of the write.

[1.0.0rc1] - 2026-07-31

First public release candidate. Magνs was developed privately up to this point, so everything below is new to anyone outside the project. The entries are still grouped as Added/Changed/Fixed/Removed, describing the development history that produced this release – “Changed” and “Fixed” are relative to earlier private states of the code, not to any published version – because that history is the most useful record of why the code looks the way it does.

Added

• Notebook 09 is now a Lorentz-invariance-violation notebook. It had been a scratch pad: no markdown, no figures, and no LIV content at all – just print(module.__all__) and a print(sys.path) that had stopped working. It is now built around the one thing that makes LIV findable, the energy scaling: the vacuum term falls as 1/E and the matter term is flat, so both switch off at high energy, while the LIV term grows as E^n and switches on. Four figures follow from that – where standard oscillations stop and LIV keeps going, how the operator dimension n_liv sets that crossover, that matter does not rescue the standard prediction (V_CC does not grow with energy either), and how a null result at high energy becomes a limit on the LIV eigenvalue, with the b*L*E^n = 1 estimate landing where the curve visibly departs.

• A magnus.plotting module of pre-packaged figures, so that a plot in the notebooks costs one call rather than the twenty-five to forty lines of gridspec_kw, tick locators, legend keywords and savefig that each figure used to carry. Taking stock of the fifty-odd notebook figures first showed that most are the same figure with different data: curves against baseline, against energy, against a sterile mixing angle, and the matrix-exponential convergence studies all reduce to one shape – curves against a swept variable over an optional relative-error subpanel – which is plot_curves, with plot_probability_vs_baseline and plot_probability_vs_energy as presets over it. Only three layouts are genuinely distinct and get their own functions: plot_probability_with_profile, plot_biprobability and plot_oscillogram, plus plot_probability_with_average for the decohered overlay. Defaults reproduce the existing house style exactly, and every function returns (fig, ax) so a packaged figure is a starting point rather than a dead end. prob_label absorbs the helper that had been copied into several notebooks, extended to cover the sterile states those notebooks needed but it did not.

• plotting.plot_curves_stacked, for small multiples. Auditing which notebook figures still built their own axes turned up one shape the module had missed: the same plot repeated once per case down a shared abscissa – notebook 07’s four panels, one per detector. It is worth packaging for the same reason as the rest: the reader compares between panels, so every panel must carry identical limits, scales and tick spacings, while the abscissa labels, title, legend and ordinate label belong to exactly one panel each. Hand-built that is four formatting loops plus a frameless full-figure subplot added purely to hang a shared label on, and it is where a stack quietly stops being readable once one panel drifts. legend_proxies also retires a trick the notebook used: plotting dummy points outside the axis limits to manufacture legend handles for entries that describe a line style rather than any one curve. Notebook 07’s two remaining hand-built figures now use this and plot_curves; the second turned out to need no new function at all, only a gridspec_kw that had never done anything.

• osc_prob_energy_baseline(..., cumulative=True): a whole baseline scan from one traversal of the profile. The evolution operator is a time-ordered product, so U(0->L2) = U(L1->L2) U(0->L1): every requested baseline is a prefix of the next, and recording the running product yields the whole scan at once instead of re-walking the profile N times. It is the reduce already in osc_prob with its intermediates kept rather than discarded.

  The grid is the union of the requested baselines (so each answer lands on a slab edge and is read off, never interpolated), a uniform accuracy grid, and any t_breakpoints. Sizing that accuracy grid is the one way a cumulative scan goes silently wrong — the traversal has nothing to compare itself against — so it is not guessed: one ordinary adaptive osc_prob call at the longest baseline reports the slab count it needed, which is the definition of the accuracy grid, and brings the existing safeguards and warnings with it. On a solar profile a plausible-looking guess of 2000 slabs is wrong by 1.6e-2 where the inherited number is right.

  Chunked traversal and conversion to probabilities at each snapshot are requirements rather than optimizations: they keep peak memory at O(block) + O(result) instead of holding N complex unitaries beside the answer.

  Measured against solve_ivp: a 1000-point solar scan takes 12.0 s per-point for an error of 5.6e-5, and 0.10 s cumulative for 5.1e-6 — 124x faster and 11x more accurate. Opt-in rather than automatic: the two paths use different grids, so results differ within the requested tolerance.

• Matplotlib is an optional dependency, declared as the plot extra (pip install 'magnuspy[plot]'). The engine still needs only NumPy, SciPy and joblib: someone computing probabilities inside their own analysis code should not have to install a plotting stack. magnus.plotting imports cleanly without Matplotlib – it defers the import into the calls that draw – so import magnus works on a core-only install and only a plotting call raises, as MatplotlibNotFoundError, naming the command to fix it.

• The PyPI distribution is named magnuspy; the import package remains magnus. Plain magnus was already taken on PyPI by an unrelated project, so pip install magnus would have fetched someone else’s package and the release workflow would have failed on upload. The two names are independent in Python packaging, so pip install magnuspy then import magnus is all that changes, and the console script is still magnus. version.py looks the distribution up under the new name: querying the import name would have raised PackageNotFoundError and fallen through to parsing pyproject.toml, which is absent from an installed wheel, so every installed user would have reported 0.0.0+unknown.

• A license. Magνs is released under the GNU General Public License v3.0 only; LICENSE carries the full text, and it is declared in pyproject.toml as the PEP 639 SPDX expression license = "GPL-3.0-only", which is what appears in the built distribution’s metadata as License-Expression. It previously read TBD, so a pip install would have shipped the license text while declaring no license at all. The build requirement moves to setuptools>=77, the first version that understands the SPDX expression form; note that PEP 639 forbids pairing it with a License :: classifier, so there is deliberately none. Referenced from the README, the docs landing page, and both file trees.

• Status badges on the README and the docs landing page: CI tests, code quality, documentation, license, the supported Python version (3.10+, matching requires-python and the CI matrix rather than a copied-in default), the code style (ruff, which lint.yml actually enforces), and PyPI downloads via pepy.tech. The downloads badge reports on the distribution name, magnuspy, and stays blank until the first release is published – pepy serves a 404 for a project PyPI does not know, so it renders as a broken image rather than as a count of zero.

• magnus.expansionterms: derives the terms of the Magnus expansion from the Bernoulli-number recursion symbolically, at any order, in exact rational arithmetic (bernoulli, bernoulli_factor, omega_terms, magnus_terms, count_terms, format_term, print_magnus_terms). The numerical core’s coefficients were typed in and nothing checked them against the recursion they come from; the test suite now regenerates them and compares, order by order, at machine precision. See docs/source/expansion_terms.rst.

• Magnus orders 7 through 10. MAGNUS_EXP_ORDER_MAX rises from 6 to 10, and is now defined once in magnus.magnus and re-exported by globaldefs, which used to carry its own copy of the number. The default expansion order is unchanged. Orders 1-6 keep their hand-written expressions, which are hot and worth reading; 7-10 are generated from the closed form of the recursion (every term is a right-nested chain of lower-order Omega_m around A, indexed by the compositions of n-1), with shared suffixes memoized so each distinct nested commutator is built once. Verified two ways: exact agreement with the symbolic generator, and measured convergence rates that keep improving with order (order 8 reaches ~h^10 against an ODE ground truth where order 6 reaches ~h^8).

• magnus.magnus.MagnusHighOrderCostWarning, raised when an order above 6 is requested with a quadrature method. The number of terms roughly doubles per order, and the measured cost per slab is 2.7x order 6 at order 7, rising to ~17x at order 10. Higher order does converge faster in the slab width, so it is a trade rather than a mistake – but narrowing the slabs at order 4 or 6 often reaches a given accuracy for less total work.

• Command-line calculator (magnus prob, also runnable as python -m magnus) for computing a single oscillation probability from the shell, covering vacuum, matter (constant/exponential density), Earth, and Sun, with standard, NSI, and LIV scenarios, for 2-5 flavors. See docs/source/cli.rst.

• Full pytest suite (tests/) and GitHub Actions CI: tests.yml (matrix across Python 3.10-3.12), lint.yml (Ruff), pages.yml (Sphinx -> GitHub Pages), publish.yml (PyPI on release).

• Test-coverage measurement. pytest-cov joins the test extra, the settings live in [tool.coverage.run] in pyproject.toml, and a separate Coverage job in tests.yml reports the figure on each run’s summary page and uploads coverage.xml as an artifact. Branch coverage is enabled: line coverage alone overstates how well this package is tested, because oscprob.py is largely thin wrappers that one parametrized test sweeps in a single pass, whereas what matters is whether the dispatch chain, the refinement caps and the warning paths are each taken in both directions. The build fails below 90%, a floor set two points under the measured 92% so that it catches a regression without tripping on the fraction of a percent that moves between interpreters; it was left unset until the figure had settled, since a threshold invented before the first measurement either sits below the real figure and never fires, or above it and blocks unrelated work. The job also carries a Codecov upload step that stays dormant until a CODECOV_TOKEN secret exists, so nothing leaves the repository until public coverage reporting is deliberately switched on.

• Two structural sweeps closing what that first coverage run found: 23 of the 36 NSI/LIV osc_prob_* wrappers, and 25 of the hamiltonian_* builders the package exports, were executed by nothing at all – not by the tests, and not by the library either, which reaches for the *_energy_independent variants instead. Several of the wrappers even appeared in a parametrize list, but in tests that inspect a signature or the source text without ever calling the function, so a mistyped keyword in any of them would have shipped unnoticed. test_every_bsm_wrapper_runs_and_is_unitary now calls every one of them with non-zero NSI/LIV parameters and checks the result is a valid probability matrix, and test_every_exported_hamiltonian_builder_is_hermitian builds every exported Hamiltonian once, with complex off-diagonal couplings, and checks Hermiticity. Both discover their own subjects – from the module and from __all__ respectively – so a name added later is swept without anyone remembering to extend a list. Both were verified to fail against a deliberately broken library before being kept.

• Tests for the paths that report failure rather than a result: the refusals of the closed-form interaction-picture integrator, both ways adiabatic.hybrid_propagator can decline to certify a non-converged patch, and what the dispatch layer then does with an uncertified point – strategy='auto' abandoning the batch for the general Magnus path, strategy='hybrid' keeping the answer and raising HybridCertificationWarning. None of these had ever been executed, which is a poor place for a blind spot in a package whose central claim is that it knows when it cannot verify its own answer. Also magnus.version’s two resolution routes, which must agree; the CLI’s argument-error paths; the unknown-location error in earth.coordinates_of_named_location; the verbose=2 banner and run-parameter dump; and globaldefs.set_color_output.

• magnus.avgprob, and the average keyword that reaches it from every oscillation-probability function through the shared **kwargs chain: the phase-averaged (fully decohered) probability, which is the exact L/E -> infinity limit astrophysical neutrinos arrive in. For a position-independent Hamiltonian it is closed-form, one eigendecomposition rather than the resolution of some 10^15 radians of phase, and for vacuum it depends on neither energy nor baseline, so one matrix serves an entire flux calculation. A position-dependent profile decoheres in the eigenbasis at production and is carried along the levels of the instantaneous Hamiltonian, with level-crossing probabilities taken from the convergence-checked Magnus patches of magnus.adiabatic rather than from a Landau-Zener formula – they reproduce that formula to a few parts in a thousand, which is the check rather than the method. A profile with discontinuities (PREM) has no closed form and is averaged over an energy window instead, which is a different quantity and says so. Whether the averaged limit applies at all is decided per pair of eigenvalues from the phase they accumulate, not assumed: pairs that have neither decohered nor stayed coherent raise magnus.oscprob.PhaseAveragingWarning, and near-degenerate ones are summed coherently within blocks, since the naive incoherent sum returns a spurious mixture where the correct answer is that nothing oscillates. Documented in docs/source/averaged_probability.rst, with worked examples in notebooks/10_magnus_averaged_probability.ipynb.

• magnus.oscprob.IP_EXP_N_SLABS_CAP and magnus.oscprob.IP_EXP_LOOP_CAP: the interaction-picture integrator’s slab and loop ceilings, previously written as bare numbers inside the function. Naming them changes no behavior, and makes the method’s conduct at the ceiling testable at a small cap – reaching two million slabs to observe what happens at the boundary costs gigabytes and minutes, so with the values inlined those paths could not be tested at all.

• docs/check_doc_snippets.py, which executes every .. jupyter-execute:: block in the documentation – in the RST pages and in the docstrings autoapi renders – and reports the page, line and traceback of any that fails. The documentation’s examples are already run at build time, so a broken one fails CI; what this adds is the ability to find out in about a second, rather than from a full Sphinx build. The fast build used while writing docs stubs those directives out, so it validates the prose and the cross-references while saying nothing at all about the code, and a page can build cleanly while being broken.

• tests/test_validation.py, covering the guards that reject bad input: the per-flavor parameter dictionaries, flavor indices and expansion orders out of range, mismatched energy/baseline arrays, negative densities and composition ratios, and the Earth entry point’s location, zenith-angle and baseline combinations. Each asserts that the input is refused as a ValueError naming the parameter at fault, rather than as whatever the interpreter happened to raise further downstream; writing them is what uncovered the except KeyError bug listed under Fixed.

• Three branches that cannot be reached are now marked # pragma: no cover, each with the argument for why written beside it: the is_repeat return in the interaction-picture integrator (below the ceiling the slab count strictly increases, at the ceiling every branch returns, so a repeat can never be compared), its no-progress guard (round(2n) == n has no solution for n >= 1, and it exists only to make a smaller growth factor safe), and the CLI’s solar baseline check (the general one rejects that input first). The Windows-only ANSI setup line in globaldefs is marked too. They are kept rather than deleted: each becomes live again if a nearby constant changes.

• A sweep asserting that both construction routes of every builder offering the choice agree: the hardcoded expression for each matrix entry, and the explicit product of the mixing matrix, the mass matrix and the conjugate transpose that it was derived from. Nothing had been keeping the two from drifting apart for the Hamiltonian builders – the existing checks covered only the 4x4 and 5x5 mixing matrices – and the product route was unexecuted for both vacuum energy-independent Hamiltonians, which oscprob calls on every run. Antineutrinos are included, since conjugating the mixing matrix is a separate line in several builders.

• Gauss-Legendre commutator-free integrators (integration_method='gl'), silent vectorization of Hamiltonian/density-profile evaluation, an energy-batched scan engine for separable Hamiltonians, adaptive slab refinement with warm-starting across scan points, and slab edges aligned with PREM layer boundaries.

• osc_prob_earth and osc_prob_sun: generic entry points that accept an arbitrary user-supplied Hamiltonian in the Earth/Sun environments.

• Sphinx documentation (published to GitHub Pages) covering installation, quick start, the CLI, the full function listing, code architecture, methodology, and this changelog.

• Root-cause regression tests for every bug listed under Fixed below.

• magnus.globaldefs.NUFIT_GLOBAL_FITS: best-fit standard three-flavor oscillation parameters from every NuFit global-analysis release, v1.0 (2012) through v6.1 (2025), by mass ordering and by each release’s secondary category (with_SK/without_SK for v4.0+, LEM/LID for v2.1, free_fluxes_rsbl/huber_fluxes_no_rsbl for v1.0-v1.3). Transcribed directly from the official parameter tables at nu-fit.org.

• magnus.globaldefs.load_nufit_params(version, ordering, category): loads a specific release/ordering/category as a plain {s12, s23, s13, dCP, D21, D31} dict, directly usable as keyword arguments to any osc_prob_3nu_* function.

• A “When is Magνs not the right tool?” section in the README and docs, covering quantum decoherence, open-system/bath coupling, neutrino decay, and self-consistent collective oscillations.

• A short pronunciation note at the top of the README and of the docs landing page: Magνs is said just like “Magnus”, with the neutrino symbol ν standing in for the “nu” syllable (with a nod to the Danish pronunciation, given where most of it was written).

• magnus.matter.exp_density_profile(density_matter_central, l_scale): builds an exponential density-profile callable tagged for the new fast interaction-picture integrator (see below); osc_prob_2nu_matter_exp_density and osc_prob_2nu_sun (and their NSI counterparts) now build their profile through it instead of an untagged lambda, so they pick up the speed-up automatically. A fast, closed-form interaction-picture Magnus integrator (magnus.oscprob._osc_prob_ip_exp_dispatch/_osc_prob_ip_exp_core) for two-flavor oscillations in a genuine exponential matter profile (the Sun): it factors the (possibly huge, at low energy) constant vacuum phase out of the Magnus expansion analytically, in closed form, instead of resolving it slab by slab, leaving only the matter-potential envelope – whose exact exponential integral is also closed-form – to be integrated. This directly fixes the bottleneck described under “Fixed” below for the two-flavor Sun wrappers: osc_prob_2nu_sun/osc_prob_2nu_matter_exp_density (and their NSI counterparts) now return correct, warning-free probabilities in a fraction of a second across the realistic solar-neutrino energy range (~0.1-18 MeV), where the general slab-refinement method previously either took many seconds to fail or returned a silently inaccurate result. The fast path is unconditionally exact in the accumulated vacuum phase and in the density profile’s shape (no local-constant or local-linear approximation); its only approximation is truncating the resulting interaction-picture Magnus series at first order, which is validated against solve_ivp across the realistic energy range and is checked at runtime (both via the size of the neglected term and via successive- refinement agreement) before being trusted, with a transparent, lossless fallback to the general method whenever it is not (e.g., near an MSW resonance, where the matter term is no longer a small perturbation on the vacuum splitting). Three- to five-flavor scenarios are not yet covered: the neglected term’s coefficient grows by roughly three orders of magnitude going from two to three flavors (one off-diagonal mass-pair vs. three, each with an O(1) diagonal mixing contribution), which pushes the slab count needed for a certified answer far past what stays fast; those wrappers therefore still use the general method unconditionally, exactly as before – a natural target for a follow-up (e.g., a genuine adiabatic/ WKB treatment, which does not have this scaling problem).

• magnus.adiabatic: a new module implementing exactly the adiabatic/WKB follow-up flagged above, generalized to any number of flavors and any Hamiltonian (not just the 2-flavor case the interaction-picture fast path above is restricted to). adiabatic_propagator(H_func, l0, l1) computes the evolution operator via pure adiabatic (instantaneous-eigenbasis) transport – a dynamical phase (Simpson-integrated) plus a geometric (Berry) phase, captured implicitly via discrete parallel transport of the eigenvectors, with no restriction on the accumulated phase or on the Hamiltonian being real. find_resonance_candidates/ find_nonadiabatic_windows(H_func, l0, l1, threshold=0.1) locate every position where adiabatic transport breaks down, for any pair of levels, via an exact Hellmann-Feynman diagnostic (no eigenvector finite differencing, which is gauge-ambiguous): candidates are exact critical points of each pairwise eigenvalue gap, and the adiabaticity verdict is the Landau-Zener-like gamma_jk = |<v_j|dH/dl|v_k>| / gap_jk^2; windows are grown to their physical width (independent of the search-grid spacing) and merged when they overlap, so any number of simultaneous or sequential resonances are handled uniformly. hybrid_propagator(H_func, l0, l1, rtol=1e-3, atol=1e-3) composes adiabatic transport with an exact local Magnus patch (magnus.magnus.magnus_expansion_multislab) at each non-adiabatic window, stitched via the exact composition law of quantum evolution (so the result is exactly unitary regardless of approximation accuracy), and self-certifies by tightening the adiabaticity threshold and two internal grid densities together until two successive results agree. Validated against solve_ivp on standard and BSM (NSI) 3-, 4-, and 5-flavor Hamiltonians, and on synthetic cases with two independent and two merging resonances: 0-window cases match to ~2e-4 at 3,600x-4,800x the speed of solve_ivp; 1-2-window (patched) cases match to 9e-4-2.9e-3 at 30x-91x the speed. See docs/source/adiabatic_strategy.rst for the full derivation and validation.

• strategy parameter ('auto' (default), 'hybrid', or 'magnus') on osc_prob_matter_std_potential, osc_prob_matter_nsi, osc_prob_liv, and every wrapper built on them (every osc_prob_{2,3,4,5}nu_sun, osc_prob_{2,3,4,5}nu_sun_nsi, and osc_prob_{2,3,4,5}nu_sun_liv function explicitly; every other wrapper via **kwargs), and also on osc_prob_sun/osc_prob_earth (the fully generic, arbitrary-Hamiltonian entry points, via a new _osc_prob_hybrid_dispatch_generic/ _osc_prob_with_potential code path, since these bypass osc_prob_matter_std_potential/_nsi/osc_prob_liv entirely). 'magnus' reproduces the exact behavior from before the adiabatic strategy was added; 'hybrid' additionally tries magnus.adiabatic.hybrid_propagator for any position-dependent, breakpoint-free potential with a requested tolerance, returning a best-effort result plus the new HybridCertificationWarning if it fails to self-certify; 'auto' tries 'hybrid' first and falls back silently to the 'magnus' strategies otherwise. Unlike the interaction-picture fast path, this applies to any number of flavors and does not require a tagged exponential profile. For osc_prob_earth, 'hybrid'/'auto' almost always fall back to 'magnus' in practice, since the PREM density profile’s layer-boundary breakpoints are essentially always non-empty for a real trajectory – the general breakpoint-free requirement above, not a special case.

• magnus.oscprob.HybridCertificationWarning (subclasses ToleranceNotAchievedWarning): raised only when strategy='hybrid' is explicitly requested and the hybrid propagator fails to self-certify for at least one requested point.

• tests/test_adiabatic.py: unit tests for magnus.adiabatic (unitarity, resonance-candidate detection, window growth/merging in both directions, and solve_ivp cross-checks for 3-5 flavor Hamiltonians), plus new regression tests in tests/test_oscprob.py: test_sun_2nu_default_strategy_avoids_tolerance_cap demonstrates the fix directly through osc_prob_2nu_sun (the same (energy, baseline) point that test_tolerance_cap_warns shows still hits the refinement caps under strategy='magnus' is resolved, warning-free and matching solve_ivp, under the new default); test_generic_osc_prob_sun_hybrid_strategy_resolves_hard_case confirms the same fix through the fully generic osc_prob_sun entry point (a separate code path, _osc_prob_hybrid_dispatch_generic, not exercised by the other test); test_generic_osc_prob_earth_strategy_falls_back_to_magnus confirms osc_prob_earth is unaffected (PREM breakpoints disable the hybrid dispatch).

• notebooks/12_magnus_adiabatic_hybrid_strategy.ipynb: live comparison of strategy='auto'/'hybrid'/'magnus' for 2- through 5-flavor Hamiltonians (standard oscillations and an engineered BSM/NSI resonance), each cross-checked against a tight-tolerance solve_ivp ground truth in both runtime and accuracy, plus a real-data plot of the instantaneous eigenvalues and detected non-adiabatic window for the 3ν BSM case. Reproduces, live, the validation described in docs/source/adiabatic_strategy.rst.

Changed

• integration_method='gl' now raises for orders above 6 rather than silently computing an order-6 result. The Gauss-Legendre commutator-free schemes are separately derived integrators, not products of the Magnus recursion, so they do not extend along with it. The check sits in _gl_nodes as well as in the input validator, since the validator is skipped when validate_input=False and that flag would otherwise reopen the silent-degradation path.

• max_n_slabs is now method-aware: it defaults to None, meaning “use the cap that suits integration_method” – 20000 for 'gl', 2000 for the cumulative-quadrature methods (magnus.oscprob.MAX_N_SLABS_DEFAULT). An explicit value is always used as given, so this changes nothing for callers who set it. A single cap could not serve both families: 'gl' costs 1-3 Hamiltonian evaluations per slab against the quadrature methods’ n_tpts_per_slab, so at a shared cap of 2000 it hit the ceiling on problems it could resolve comfortably – eV-scale sterile splittings over an Earth-crossing baseline need about 8,600 slabs – and emitted ToleranceNotAchievedWarning on answers that were roughly 1,600x more accurate than the quadrature methods achieved within that same cap. Even at 20000, 'gl' remains the cheaper worst case (40,000-60,000 evaluations against the ~200,000 that 2000 quadrature slabs at 100 points already permit).

• Breaking: integration_method now defaults to 'gl' (Gauss-Legendre commutator-free collocation) instead of 'trapezoid', everywhere it appears – the 13 signature defaults across magnus and oscprob, and magnus prob --integration-method. For a Hamiltonian that is smooth within each slab, which layer-aligned slab edges make the common case, 'gl' is simultaneously the faster and the more accurate choice: it needs only 1, 2, or 3 Hamiltonian evaluations per slab for orders <= 2, <= 4, <= 6, with its quadrature order matched exactly to the truncation order, where the cumulative-quadrature methods sample n_tpts_per_slab points and can let quadrature error dominate the truncation error. Numbers computed without passing integration_method explicitly will therefore change slightly, and two further behaviors switch on with it, both of which were already implemented and gated on 'gl': n_tpts_per_slab no longer participates in the adaptive refinement (accuracy is set by the slab count alone), and the physics-informed starting slab count from magnus.suggest_n_slabs is now applied by default. 'trapezoid' and 'simpson' remain fully supported and are the better choice when the Hamiltonian has a kink or a discontinuity inside a slab, where Gauss-Legendre loses its order advantage.

• Collapsed roughly 1,150 lines of duplicated refinement/logging keyword-argument declarations across ~60 wrapper functions into a single source of truth (internally called the “G1” refactor), with a permanent test (test_no_wrapper_redeclares_standard_refinement_kwargs) guarding against the pattern recurring.

• Rewrote the Magnus-expansion numerical core: corrected higher-order term coefficients, added order-6 support, and restructured the term recursion and matrix exponentiation for batched/vectorized evaluation.

• Package layout consolidated under src/magnus/ (src-layout) with proper __init__.py files; version metadata unified to a single source of truth.

• Restructured the package to be flatter: earth, globaldefs, magnus (the numerical core), matter, oscprob (the main wrapper API), and oscprobstd are now flat sibling modules directly under src/magnus/, instead of each living inside its own single-file subpackage directory. hamiltonians remains a genuine subpackage (four distinct flavor-count-specific modules). Every __init__.py (top-level and hamiltonians/) now uses explicit, named imports and a hand-written __all__ instead of from .module import * plus a dir()-computed __all__, which had been silently leaking implementation-detail names (np, Optional, Callable, …) into the public namespace.

• magnus.authors and magnus.version are no longer part of the public API surface (excluded from __all__ and from the Sphinx autoapi-generated docs), though they remain importable internally for the CLI’s --version flag and oscprob’s banner-printing.

• Docstring equations that were still written as indented plain text now render as real LaTeX. The commit that converted ~600 instances of plain-text symbol notation to :math: roles left the bare equations alone; magnus.py’s module docstring (the Magnus expansion and its Bernoulli-number recursion), MagnusConvergenceWarning’s convergence criterion, earth.prem_layer_edges_along_chord’s chord-crossing quadratic, matter.density_matter_func_exp’s profile, and oscprob._osc_prob_ip_exp_core’s two interaction-picture equations are now .. math:: blocks whose LaTeX matches docs/source/methodology.rst for the identical mathematics.

• Physical units carrying exponents are now math-mode throughout the docstrings, so they render as eV², g cm⁻³ rather than as the literal text eV^2, g cm^{-3}. The convention adopted, and applied uniformly across oscprob, globaldefs, matter, earth, and cli, is to math-wrap a whole unit group when — and only when — it carries an exponent ([:math:`\text{g cm}^{-3}`]), matching the pre-existing \text{cm}^{-3} usage in oscprob’s module docstring; exponent-free units ([eV], [radian]) stay plain text, since they already render correctly. docs/source/functions.rst’s lone :sup:-based unit was converted to match. This also normalizes the eV^-1/eV^{-1} spelling inconsistency.

• All 70 remaining >>> doctest-style prompts (in oscprob’s module docstring and in globaldefs.load_nufit_params) are now .. jupyter-execute:: blocks, the convention already used by the ~60 osc_prob_* wrapper docstrings. No doctest runner has ever existed in this project — no --doctest-modules, no doctest CI step, no doctest_namespace fixture — so these examples were never executed or checked by anything, and had silently drifted out of sync with the code (see “Fixed” below). As executed cells they now show genuine, always-current output and will fail the docs build if they ever break again.

• Breaking: invalid input now raises ValueError instead of printing a message and calling sys.exit(1). There were 62 such aborts across oscprob and earth; a library that terminates the interpreter cannot be recovered from in a notebook, a scan loop, or a caller that wants to fall back, and it made the failures impossible to assert on in tests. Most were already raise ValueError(...) inside a try block whose except immediately swallowed it and exited, so the messages are unchanged — they now propagate instead of being printed. validate_input_battery follows suit: it raises rather than returning 1, so its return type is now None (it previously returned 0/1, and every caller compared against 1). Error messages raised as exceptions use the plain-text gd.ERROR_MSG_NO_COLOR prefix, since ANSI codes are meant for a terminal and end up in tracebacks and logs. The magnus CLI catches these and reports them as ordinary argument errors, so its behavior is unchanged.

• A single source of truth for the version number: the version field of pyproject.toml. magnus/version.py now resolves it via importlib.metadata when installed, falling back to parsing pyproject.toml when running off src/ on the path, and docs/source/conf.py imports it rather than repeating it. The eleven decorative per-module __version__ strings (which disagreed with each other — "2.0" in magnus.py, "0.10.0" in oscprob.py, "1.0" elsewhere — and which nothing ever read) are gone.

• The strategy keyword is now reachable from the command line as magnus prob --strategy {auto,hybrid,magnus}, for the environments where the Hamiltonian actually depends on position (sun, earth, and matter with --density-profile exp). It was previously Python-API-only, so CLI users silently got 'auto' with no way to opt out or to force it.

• globaldefs.set_color_output(enabled) turns the ANSI color in the warning/error/tolerance prefixes on or off. WARNING_MSG_NO_COLOR and its siblings had existed and been exported since the beginning but were never used by anything: every call site hardcoded the colored variant, so there was no supported way to get clean output into a log file or a rendered notebook.

• The documentation build now runs on pull requests, not only on pushes to main. It was previously built only by pages.yml, so a page that failed to build, or silently dropped out of the sidebar, was not caught until after it had been merged. The pull-request build is deliberately stricter than the deploy build – -W turns Sphinx warnings into errors, which is what catches a broken cross-reference or a page missing from the toctree – while the deploy build stays permissive, so a late warning can never be what stops a release from publishing.

• ruff check is now blocking in CI rather than continue-on-error, with the rule configuration in [tool.ruff.lint] in pyproject.toml, including an explicit select. Leaving the selection to ruff’s default is not safe for a blocking check: the default widened in ruff 0.16 and turned CI red on a codebase that had not changed. Two codebase-wide conventions are exempted explicitly (E741, since l is the standard symbol for position here, and E701 for the one-line cleanup guards); everything else was fixed, so the tree is clean and a new finding fails the build instead of being reported into a green checkmark.

Fixed

• A batched solar call could exhaust the machine’s memory. _osc_prob_ip_exp_core — the closed-form interaction-picture integrator behind osc_prob_2nu_sun and its NSI/LIV siblings — built temporaries of shape (n_energies, n_slabs, d, d) while its ladder doubled n_slabs toward IP_EXP_N_SLABS_CAP = 2_000_000. The reasoning recorded beside that ceiling was about time only (“each slab costs one 2x2 eigendecomposition”), and never accounted for the energy count multiplying the working set. Measured, on the documented solar use case: ~1.3 GB per energy — 1.56 GB at one energy, 5.34 GB at four, and a MemoryError of shape (8, 2000000, 2, 2) at eight. Nothing about the call is pathological; it is the advertised batched form of a public wrapper, and notebook 03 scans 1000 energies over exactly that range. The notebooks never hit it only because they use raw osc_prob loops. The working set is now tiled over both axes against a fixed budget (BATCH_WORKING_ENTRIES), giving a peak that is flat in the energy count — 79 MiB at 4 energies and 79 MiB at 256. The tiling is exact: the slab product is folded in the same order with the same parenthesis nesting, and a test pins the output at bit equality against an untiled run. Found while measuring the energy axis, and pre-existing — reproduced unchanged at 155e01e.

• The same integrator burned twenty-one refinement levels to reach refusals it could have predicted. Certification requires max|Omega_t| below a trust threshold, and that maximum is bounded below by the diagonal entries, which have a closed form. When even that bound exceeds the threshold at the slab ceiling, no reachable slab count can certify. Two evaluations of the potential and no allocation now detect it, and the method refuses immediately instead of doubling its way to the same answer. It is a bound rather than an estimate, so it can only report “impossible” and cannot abandon a case that would have converged. At 10 MeV over a solar radius it fires for any tolerance of 1e-4 or tighter, or any baseline beyond two solar radii.

• A scan whose result could not fit reported it as an out-of-memory kill. Tiling bounds the engines’ working set, but nothing shrinks the answer: N points over d flavors is N*d*d floats either way. osc_prob_energy_baseline now checks that against the operating system’s free-memory figure and raises a MemoryError naming the size, rather than letting an overcommitting kernel take the machine down instead of the process. The check costs one multiply below a 64 MiB floor, and never blocks where free memory cannot be read.

• A requested n_slabs was silently discarded whenever a tolerance was on, and that certified wrong answers. osc_prob documented and implemented “if rtol or atol is given, n_slabs is ignored”: the adaptive ladder started at min_n_slabs = 1 regardless of what the caller asked for. The seed that replaced it, magnus.suggest_n_slabs, measures the integral of the Hamiltonian along the path, which is blind to structure that averages out. On notebook 03’s castle-wall profile – 50 square density walls – the whole trajectory accumulates only ~9 radians, so a call asking for n_slabs=150 was seeded with 2 slabs and stopped at 4. Four slabs cannot see fifty walls, and the ladder they sit on does not converge, it thrashes: 0.43, 0.13, 0.13, 0.64, 0.12 at 2, 3, 4, 5, 6 slabs. The successive-iterate test fired on the accidental 3-vs-4 agreement and returned a probability wrong by 0.855, with no warning. Tightening rtol does not help – the comparison is between two answers that both failed to see the profile. n_slabs is now a floor on the ladder, in osc_prob and in both batched scan engines: refinement starts at max(min_n_slabs, n_slabs) and only ever climbs, clipped at max_n_slabs so a floor above the cap raises the existing not-achieved warning instead of stepping the ladder back down. With the default n_slabs=1 the floor is inactive and nothing changes. The regression test is built on a solve_ivp oracle, not on agreement between two osc_prob grids – that kind of agreement is what let this through.

• Notebooks 02 and 03 shipped castle-wall figures drawn from those wrong probabilities. Across a 6000-point baseline scan, 7.2% of points were off by more than 1e-2. The scans now get their n_slabs=150 honored, and they also pass the wall positions as t_breakpoints: the profile is a step function, and high-order quadrature reaches its nominal order only when the Hamiltonian is smooth inside each slab. Against a converged reference the worst point improves from 0.855 to 1.0e-3 – and the scans run faster than the wrong version did, because slab edges placed on the discontinuities buy more accuracy per slab than piling on uniform slabs.

• Eight cells in notebooks 02 and 03 taught a performance lesson that is not true. They asserted that an osc_prob call costs “~4.5 ms, almost entirely fixed entry-path cost”, evidenced by “n_slabs=1 and n_slabs=150 both measure 4.5 ms per call, which is how you can tell the physics is not what costs”, and predicted that looping instead of passing an array would turn a 2 s cell into 45 s. Measured: 0.26-0.46 ms per call, and on a profile that actually varies the same call costs 0.34 ms at n_slabs=1 against 10.4 ms at n_slabs=2000 – the integration is nearly all of it. The loop-vs-array ratios are 1.4x-2.3x, not the ~22x claimed. Passing the array is still the right advice and the cells still do it; the numbers and the reasoning behind them are now the measured ones.

• Every notebook that used matter effects was computing vacuum. Notebooks 01-10 built the coherent forward potential by calling matter.num_density_e_func with a density in g cm^-3 but without density_matter_is_in_g_per_cm3=True – for constant densities and for earth.density_matter_func_prem alike. That yields VCC = 8.8e-32 eV instead of 3.8e-13 eV, a factor of 4.3e18, so the matter term was ~20 orders of magnitude below the vacuum one. In notebook 02 the “matter” probability came out bit-identical to the vacuum one. 19 call sites.

• The matter potential was subtracted where it should be added. For neutrinos the library computes H_vac + h_matt(VCC); the antineutrino sign flip lives inside VCC (matter.vcc_func_from_rho_func applies s = 1 if not nubar else -1), and matter.VCC_func – which the notebooks used – always returns a positive potential. 45 sites in notebooks 02, 03, 04, 07 and 08 carried a leading minus on a neutrino Hamiltonian. Checked against the closed-form oscprobstd result: the corrected sign agrees to 2.4e-14, the old one is wrong by up to 135%. These two bugs masked each other – with the matter term 20 orders down the sign was invisible, and the notebooks’ own relative-error subpanels read a healthy 1e-12 because the standard formula was being fed the same wrong potential.

• PREM was sampled at the center of the Earth. Notebooks 04, 05 and 07 wrote VCC_func_prem(r/gd.CONV_KM_TO_INV_EV) where r was already in km, evaluating the profile at r ~ 1e-6 km and so using a constant 13.09 g cm^-3 everywhere instead of the layered profile. The same notebooks used the correct VCC_func_prem(r) for their density panels, so the plotted profile disagreed with the physics behind it. Cross-checked against osc_prob_3nu_earth: the corrected Hamiltonian agrees to 3.2e-6, while the old one gave P(nu_mu -> nu_e) = 0.0001 against a true 0.0047. 7 sites.

• A matter term discarded as dead code. Eight Hamiltonians across notebooks 05 and 07 read return H_vac(...) followed by a bare + H_matt on the next line, which Python evaluates and throws away, so those Hamiltonians were pure vacuum. Joined onto the return.

• Notebook 07 had never run top to bottom: cell 14 referenced baseline, which is first assigned three cells later. Five figure titles in notebooks 07 and 08 named the swept variable instead of the fixed one.

• A matter density that has already been converted to natural units, but is then declared to be in g cm^-3, is now flagged with matter.DensityUnitWarning instead of being converted a second time in silence. The two scales do not overlap – the densest matter anyone models is some 1e15 g cm^-3, while any density from water upwards becomes 4.3e18 or more once converted – so the check has three orders of magnitude of margin. It is worth having because the consequences do not look like a unit error: the matter term swamps every other scale, nu_e becomes an exact eigenstate, and the calculation returns a perfectly self-consistent P_ee = 1, which reads as a broken formula rather than as bad input.

• Passing a whole entry of globaldefs.OSC_PARAMS_PREDEFINED to a probability function, as **OSC_PARAMS_PREDEFINED['OSC_PARAMS_DEFAULT'], now raises a ValueError naming the two offending keys and pointing at globaldefs.load_nufit_params. Those entries carry name and description strings alongside the six mixing parameters; unchecked, they traveled the shared **kwargs chain until magnus_expansion_multislab rejected them, naming the one function in the chain with nothing to do with the mistake. The check sits in the four middle-layer functions rather than further down, because the averaged and ordinary paths diverge before the Magnus core is reached – a guard placed later caught the strings on one path and ignored them on the other.

• matter.vcc_func_from_rho_func raised TypeError: 'float' object is not callable when given a constant electron number density together with density_is_of_number_of_electrons=True. VCC_func evaluates whatever it is handed at a position, and the sibling branch for a matter density wraps its constant in a function for exactly that reason; this one passed the bare number through. The documented way to supply a fixed electron number density therefore did not work at all. Found by the first test to reach the line.

• The two guards in validate_input_osc_prob_earth that reject a malformed loc_ini/loc_fin caught KeyError, which unpacking a sequence never raises, so neither could fire. A three-entry coordinate escaped as too many values to unpack (expected 2) and a non-iterable as TypeError – the latter breaking the convention, established across the rest of the package, that bad input raises ValueError with a message naming the parameter at fault. They now catch TypeError and ValueError. Found by writing the first test that ever reached them.

• osc_prob raised UnboundLocalError instead of returning a probability when max_num_loops < 1 was passed together with validate_input=False (the validator rejects that combination otherwise). The refinement-limit checks at the top of the loop could return P before the first iteration had produced one; they are refinement limits and only mean anything after a loop has run, so they are now guarded on that. The dead iterate_over_magnus_exp_order dispatch had been assigning P early, which hid the problem from static analysis until it was removed.

• Two mixing-matrix formula bugs (mixing_matrix_4x4 and mixing_matrix_5x5) that invalidated every sterile-neutrino (3+1, 3+2) calculation.

• hamiltonian_2nu_nsi’s eps_aa parameter was a silent no-op: it sat on both diagonal entries, making it a pure multiple of the identity (an unobservable global phase) with zero effect on any probability.

• osc_prob_5nu_matter_nsi_exp_density called the non-NSI inner function instead of the NSI one.

• A sign error in hamiltonian_2nu_liv_energy_independent’s off-diagonal term.

• Several _nsi_td/_liv_td position-dependent convenience functions crashed with TypeError on any call.

• unpack_nsi_params_from_dict/unpack_liv_params_from_dict silently returned None instead of raising for unsupported flavor counts.

• Missing nubar parameter on several matter/LIV wrapper functions.

• A boundary bug incorrectly rejecting rho_central == 0.0 in the NSI exponential-density wrappers.

• hamiltonian_3nu_liv crashed with TypeError on every call (it forwarded an incomplete argument list to its own energy-independent helper).

• Dead/unreachable code and stale, copy-pasted docstrings (wrong flavor count, wrong matrix dimensions, description of the wrong scenario) across hamiltonians{2,3,4,5}nu.py and oscprob.py, found while writing complete docstrings for every function.

• NumPy 2.0 compatibility (removed deprecated type aliases).

• ~47 docstring/type-annotation mismatches across the codebase, found by a systematic audit comparing every function’s actual return/parameter type against what its docstring documented (not just the numpy.ndarray-documented-as-list case originally flagged in hamiltonians{2,3,4,5}nu.py and oscprobstd.py, but every such disagreement anywhere in src/magnus). Also fixes two unrelated copy-paste errors caught along the way in oscprobstd.J()’s docstring (a mislabeled parameter description and a wrong worked example).

• Three latent errors in oscprob’s module-docstring examples, exposed by converting them to executed cells: np.array([1.0, 10.0 100.0]) was missing a comma (a SyntaxError), osc_prob_3nu_matter_constant_density was called without its oscprob. prefix (a NameError), and the flavor indices were listed as “NUE, NUMU, and NUMU” instead of NUTAU.

• adiabatic.hybrid_propagator could report certified=True without having certified anything. Its three refinement knobs saturate at different iterations (n_probe at 5, n_points at 6, threshold at 11), so by iteration 12 all three were pinned: that iteration recomputed bit-identical inputs and the agreement test compared a result with itself. The loop now stops as soon as every knob has saturated and reports certified=False, which is the honest answer. Covered by a regression test that fails against the old code.

• adiabatic.find_resonance_candidates (and everything built on it) no longer evaluates the user’s H_func outside the requested [l0, l1]. The finite-difference stencil reached to l0 - h and l1 + h at the endpoints, which can raise or return nonsense for a Hamiltonian defined only on its physical domain – earth.density_matter_func_prem, for one, raises beyond EARTH_RADIUS. The stencil is now one-sided at the boundaries.

• The speedup chart in docs/source/adiabatic_strategy.rst disagreed with the validation table directly above it: it showed a 25,800x bar for the case the table reports as ~30x (25,800x was a different, unlisted measurement). Both now come from VALIDATION_GRID in the new docs/make_figures.py, so they cannot drift apart again.

• docs/source/installation.rst claimed src/magnus/ had to be on the Python path “for a few modules that resolve sibling imports directly”. No such import has existed since the package was flattened; only src/ is needed. The README and quickstart.rst also still carried sys.path.extend(['src', 'src/magnus'])  # until pip packaging lands, long after packaging landed.

• HybridCertificationWarning is exported from magnus.oscprob.__all__. It was the only public class or function in the package missing from its module’s __all__, despite being raised, documented, and cross-referenced.

• Four Returns sections said np.narray instead of np.ndarray, and earth.coordinates_of_named_location’s message said “the given name of the the location”. Several other error messages had typos (wil, lengh, “only of the two”, and a mangled function name in values_to_unspecified_osc_params).

• Every :func:/:class: cross-reference in the docs now resolves. Bare names in quickstart.rst and adiabatic_strategy.rst are qualified; sibling-module references in hamiltonians{2,3}nu.py are qualified; and references to private helpers (for which autoapi never emits targets, so they rendered as plain text) are now inline literals, which is what they should have been.

• Undocumented parameters: tol in earth.density_matter_func_prem and earth.earth_radial_distance_from_depth, and A_eval_mode in magnus.magnus_expansion and magnus.magnus_expansion_multislab.

• osc_prob and the four osc_prob_{2,3,4,5}nu_vacuum wrappers documented their parameters without types, unlike the other 167 functions in the package, so the rendered docs dropped the type column for the single most important function in the API. All 99 entries now carry the type from the signature.

• Stale “Routine listings”: globaldefs was missing load_nufit_params and matter was missing exp_density_profile. cli.build_parser and cli.main were the only public functions with no .. versionadded:: tag.

• .github/workflows/tests.yml only ran on pushes to main and dev, so every push to dev-plotting – the branch this work happens on – got no CI at all. It now covers dev-* too, and installs the package with pip install -e . rather than only its requirements, so the console script and the version lookup are exercised the way a user’s install is.

• docs/source/cli.rst quoted magnus prob --help from a hand-copied paste that had gone stale (it never mentioned --version). It is now generated by docs/regen_cli_help.py, with a CI job that fails if the page and the parser disagree.

Removed

• docs/source/sandbox/, an untracked, unused pydata-theme experiment directory.

• 1,257 lines of ad-hoc scratch code in the if __name__ == "__main__": blocks of oscprob (1,218 lines, ~8% of the file), earth, matter, and magnus. The first three could not run at all: executing those modules as scripts fails at import under the flattened package layout. magnus’s block did run, but it printed an unasserted comparison of the three integration methods that tests/test_magnus_expansion.py already parametrizes over with real assertions. The only __main__ blocks left are the two genuine console entry points, cli.py and __main__.py.

• All .. versionchanged:: directives. This being the first public release, there is no earlier published version for behavior to have changed from, so the six of them described how the code works rather than what changed. Their content is kept, reworded as .. note:: blocks on the same functions.

• The eleven decorative per-module __version__ strings, and the separate [Unreleased]/[0.10.0] changelog sections, which are consolidated here: nothing was published before this release, so the split served no reader.