Magνs: Neutrino Oscillations via the Magnus Expansion

Hint

How do I say that? Just like the name Magnus — the Greek letter ν (nu), the neutrino’s symbol, simply stands in for the “nu” syllable. (And since most of this package was written while the author was based in Denmark, you are equally welcome to say it the Danish way.)

Important

Important Links:

What it can compute, with code
GitHub Repository
Example Notebooks (see also Tutorial notebooks for a guided tour)
How to cite
Changelog

Magνs computes neutrino oscillation probabilities between an arbitrary number of flavors, for any given Hamiltonian, time-dependent or -independent. Internally, it propagates the neutrino evolution operator using the Magnus expansion: rather than integrating the Schrödinger equation step by step, it exponentiates truncated time-ordered integrals of the Hamiltonian over a chain of position slabs. Any truncation of the Magnus series lives in the Lie algebra, so the resulting evolution operator is exactly unitary by construction — probabilities are non-negative and sum to one at machine precision, at any accuracy setting.

What “accurate” means here

Magνs is a numerical integrator, so unlike a closed-form method it has an error that depends on how finely it discretises. Two properties are exact regardless, and the rest is measured rather than asserted.

Exact at any order, by construction. Every truncation of the Magnus series is anti-Hermitian, so its exponential is exactly unitary — not unitary to within a tolerance. Truncating early costs accuracy, never norm.

Property	Measured agreement
Unitarity, \(U^\dagger U - \mathbb{1}\)	1e-16 to 1e-13
Magnus terms vs an independently coded Bernoulli recursion, orders 1–6	machine precision
Gauss–Legendre convergence rate under slab halving, orders 2/4/6	error ratios 4 / 16 / 64
2ν and 3ν vacuum vs the closed-form expression	machine precision
2ν constant-density matter vs the closed form, ν and ν̄	machine precision
Earth crossing (PREM) at the default `rtol = atol = 1e-3`, against a 1e-7-tolerance reference	~5e-4
Asymmetric profiles with complex Hamiltonians vs `solve_ivp`/DOP853 at `rtol=1e-12`	1e-4 to 1e-7
Energy-batched scan vs the per-point path	exactly 0.0
`n_jobs > 1` vs serial	exactly 0.0

The last two rows are the ones worth reading twice: they are bit-identity assertions, not tolerances, so an optimisation that changed an answer would fail them rather than pass quietly.

And the honest caveat. rtol/atol are a stopping criterion — the ladder halts when two successive refinement levels agree — not a bound on the error of what is returned. Usually that is conservative. It is not always: What rtol and atol actually control gives the measured detail, including a case where two levels agreed coincidentally and the answer was wrong by 0.855. Magνs warns loudly in that regime, and Implementation Details reports the measured false-alarm rate of each warning.

When is Magνs a win?

Compared to solving the propagation ODE directly (e.g., with an adaptive Runge–Kutta solver), Magνs wins when one or more of these apply:

The matter profile varies slowly compared to the oscillation length. A Magnus slab is exact for a constant Hamiltonian no matter how many oscillation cycles it spans, so the slab size is set by how fast the profile changes, not by how fast the phase winds. An ODE solver must resolve every oscillation. For a 1 GeV neutrino crossing the Earth (PREM profile), Magνs needs ~10 slabs plus the ~16 layer crossings, versus thousands of right-hand-side evaluations for solve_ivp — measured: ~2 ms vs ~360-700 ms per probability at comparable accuracy.
You scan over energy and/or direction (spectra, oscillograms, sensitivity studies). The Magnus kernel is built from fixed, data-independent matrix operations, so slabs — and, for the standard/NSI/LIV Hamiltonians, the entire energy axis — evaluate as batched NumPy/BLAS calls. Adaptive ODE integration is inherently sequential and cannot share steps across energies. Measured: a 200-energy Earth-crossing scan takes 76 ms (0.4 ms per energy); a 100×100 oscillogram takes ~2 s.
Unitarity matters more than raw local error — long baselines, small probabilities, CP/T asymmetries. Runge–Kutta iterates drift off the unitary manifold (probability leaks of ~1e-6 at typical tolerances, growing with baseline); the Magnus route has no leakage to leak, ever (probability rows sum to 1 to ~1e-14).
You want arbitrary physics with no per-model work: any number of flavors, any Hermitian Hamiltonian — sterile neutrinos, non-standard interactions, Lorentz-invariance violation, or your own matrix function of energy and position.

When is it not the best tool? For a single probability at a single energy, any method is fast enough. For extreme accumulated phases — e.g., ~10 MeV neutrinos crossing most of the Sun (~1e4 rad of matter-dominated phase) — the plain Magnus slab-refinement method can need a very large slab count, and warns (ToleranceNotAchievedWarning) instead of failing silently if it hits its caps first. This regime is now handled automatically by strategy='auto' (the default for osc_prob_matter_std_potential, osc_prob_matter_nsi, osc_prob_liv, and every wrapper built on them, including every osc_prob_*_sun* function): an adiabatic-transport-plus-Magnus-patch strategy that stays exactly unitary and is 50-25,000x faster than direct integration across the validation grid — see Adiabatic + Magnus Hybrid Strategy for the full derivation. And a tight-tolerance ODE solver remains the best reference for validation regardless — Magνs’s own test suite uses scipy.integrate.solve_ivp at rtol=1e-12 as ground truth. See Methodology for the full numerical story, including how these numbers were measured.

When is Magνs not the right tool?

Magνs solves the unitary Schrödinger equation for a Hermitian Hamiltonian: any truncation of the Magnus series lives in the Lie algebra, so the package is architecturally committed to norm-preserving, reversible evolution. That rules out several classes of problems that show up in neutrino phenomenology:

Quantum decoherence. Wave-packet separation, quantum-gravity-induced decoherence, or any model where coherence between mass eigenstates is damped over the baseline requires evolving a density matrix under a non-unitary master equation (e.g., Lindblad/GKSL), not a state vector under a Hamiltonian. Magνs has no dissipative term and cannot represent one.
Open-system coupling to a bath. Any scenario where the neutrino exchanges energy or phase information with an environment – collisional decoherence, thermal baths, stochastic scattering beyond the mean-field matter potential – needs a reduced density matrix with dissipators, which is again outside a Hermitian-Hamiltonian, pure-state framework.
Neutrino decay. Invisible or visible decay into lighter states removes probability from the system, so the evolution is no longer norm-preserving. A Hermitian effective Hamiltonian cannot encode a decay width – that requires an anti-Hermitian term, which breaks the unitarity the whole method relies on.
Self-consistent collective oscillations. Dense-environment (e.g., supernova) neutrino self-interactions, where the effective Hamiltonian depends on the (unknown, evolving) neutrino/antineutrino flavor content itself, are a nonlinear, self-consistent problem. Magνs assumes the Hamiltonian is a known function of energy and position supplied by the caller, not a functional of the solution.

If your problem needs any of the above, look instead at packages built around density-matrix/Lindblad evolution (for decoherence or decay) or dedicated collective-oscillation codes (for self-interaction problems).

What it is not

Separately from the physics above, and worth saying plainly so that nobody evaluates Magνs for a job it was never meant to do:

Not a solver for constant Hamiltonians in a hurry. It will do them, but a closed form beats an integrator every time; see When to use NuOscProbExact instead.
Not a flux, cross-section or detector code. It computes oscillation probabilities and stops there.
Not a fitting framework. There is no likelihood machinery; the probabilities are meant to be handed to whatever does that.
Not an event generator, and not an unfolding tool.

When to use NuOscProbExact instead

Magνs integrates a Hamiltonian across each slab, which is what lets it follow a density that changes as the neutrino moves. That machinery is wasted — and slower than the alternative — when the Hamiltonian does not change at all.

Reach for NuOscProbExact instead when the Hamiltonian is constant, or piecewise constant. It expands the Hamiltonian and the evolution operator in the SU(2), SU(3) and SU(4) bases, which gives a closed form rather than a numerical integration: exact up to floating-point round-off, and with no slab count to choose.

Situation	Use this	Because
Constant density	NuOscProbExact	One closed form, no integration
Piecewise constant, tens of layers — the Earth through PREM	NuOscProbExact	Each layer solved exactly, operators multiplied
Smoothly varying, slow against the oscillation	Either	Slabbing converges quickly
Smoothly varying, fast against the oscillation — the Sun, adiabatic MSW	Magνs	Slabbing needs \(\sim 10^4\) steps per resonance crossing
A shock front, a kink, a tabulated profile	Magνs	`t_breakpoints` puts a slab edge on the discontinuity
More than four flavours	Magνs	The SU(N) expansions stop at SU(4); Magνs has no ceiling
Genuinely open systems: decay, decoherence	Neither	Needs a Lindblad solver, not a unitary one

The two packages share conventions, units and parameter defaults deliberately, so a calculation can be moved between them as a cross-check. That is worth doing: agreement between two methods with different failure modes is stronger evidence than either one’s internal convergence check.

Performance

A single 3ν Earth probability takes about 2 ms at the default tolerance, and the median call across 164 Earth and solar configurations is 2 ms with the slowest at 0.90 s. Scans are what the code mostly does, and three things make them much faster without changing any answer.

Pass arrays instead of looping. Every wrapper takes an array of energies, of baselines, or both. For a position-dependent Hamiltonian the matter profile is then built once for the whole scan rather than once per point, which is what the energy-batched engine exists to do.

Write your ``H_func`` so it accepts an array of positions. The single largest factor under a caller’s control: measured at 4.6x on a 3ν exponential-density profile, with bit-identical output. A scalar-only Hamiltonian raises ScalarHamiltonianWarning once per session, naming the fix. See Writing an H_func that does not cost you a factor of five.

An Earth chord is a palindrome. A neutrino crossing a spherically symmetric Earth meets every radius twice, so the Hamiltonian is evaluated on the first half of the slab chain and the rest follows by reversal. That halves the calls to your H_func, so it is worth what your Hamiltonian costs:

Workload	Speed-up	Note
Single point, plain PREM	0.91x	a density lookup is too cheap to halve
Single point, expensive `H_func`	1.41x–1.67x
12- and 40-energy scan, expensive `H_func`	1.56x–1.64x
Energy scan, standard PREM	1.00x	the separable engine already shares the profile

magnus.magnus.USE_PALINDROME switches it off. Standard PREM scans are unaffected because the batched engine already evaluates the profile once and shares it across energies — the same saving, taken earlier.

And one cost that runs the other way. The adaptive ladder computes the probability at several slab counts and stops when two agree, so a call at a tight tolerance is doing real extra work rather than being slow. rtol=atol=None runs once at the grid you specify.

Implementation Details reports where the time goes, and what was tried and rejected.

Salient Features

Two ways to use it: as an importable Python module (the full API – see Quick Start Guide) or as a magnus command-line calculator for a single probability with no Python required (see Command-Line Calculator).
Any number of flavors, any Hamiltonian: dedicated, validated wrappers for 2ν, 3ν, 4ν (3+1 sterile), and 5ν (3+2 sterile) systems (see Available Oscillation-Probability Functions for the full listing), plus a fully generic entry point (osc_prob) that accepts an arbitrary Hermitian Hamiltonian of any dimension.
Vacuum, matter, Earth, and Sun: constant-density matter, exponentially falling density profiles, the Earth (Preliminary Reference Earth Model, including chords between named detector sites), the Sun, or any density profile you supply.
Beyond the Standard Model: non-standard neutrino interactions (NSI) and CPT-odd Lorentz-invariance violation (LIV), for every flavor count and environment above.
Magnus expansion to order 6, with the term recursion verified term-by-term against the literature, and three integration methods. The default, Gauss-Legendre commutator-free integrators, reaches orders 2/4/6 from only 1/2/3 Hamiltonian evaluations per slab; cumulative trapezoid/Simpson quadrature remains available for Hamiltonians that are not smooth within a slab.
Exact unitarity, adaptive refinement to a requested tolerance with physics-informed starting slab counts and warm starts across scans, slab edges aligned with density discontinuities, and an energy-batched scan engine for standard/NSI/LIV Hamiltonians.
Silent vectorization: Hamiltonian and density-profile functions that accept position arrays are detected and used automatically, with a safe scalar fallback.

Getting started:

Using Magnus:

How it works:

Author

Magnus was written by Mauricio Bustamante (mbustamante@gmail.com). Bug reports and questions are best raised as GitHub issues, which leave a public record others can find.

Citing

If Magnus contributed to work you are publishing, please cite it, and say which version you used – results can depend on it. How to cite has the BibTeX entry and the two or three things worth stating in the text.

License

Magνs is released under the GNU General Public License v3.0 only (GPL-3.0-only). The full text ships with the source, as LICENSE in the repository root, and inside the installed distribution.

You are free to use, study, modify, and redistribute it, including for commercial purposes, provided that derivative works are distributed under the same license and with source available. If you are unsure whether your intended use is compatible, read the license itself rather than this summary.