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.

Flexible. The Hamiltonian is an argument, not an assumption. Standard oscillations, non-standard interactions, Lorentz-invariance violation, sterile states, pseudo-Dirac pairs and a model of your own all go through the same call. Two to five flavors ship ready-made; the generic entry points take any dimension and any profile, given as a function of position.

Fast. A scan over energy or arrival direction is one batched call rather than a loop, worth one to two orders of magnitude per probability. The median call over 164 Earth and solar configurations is 2 ms; a 200-energy Earth-crossing scan takes 76 ms, and a 100x100 oscillogram about 2 s.

Accurate. Probabilities are unitary by construction at every setting, not by refinement, and on a smooth profile Magνs reaches 2.9e-13 where a composition of constant slabs floors at 2.5e-11. Where it cannot certify its own answer, it says so.

What it can compute

• Oscillations through a varying profile: the Earth’s PREM layers, any of twelve tabulated standard solar models, a supernova shock front, or any density you supply.

• The phase-averaged probability a solar or astrophysical experiment actually measures, over its energy resolution, without resolving the oscillation.

• The evolution operator itself, alongside the probabilities, for observables built from amplitudes.

• The same probabilities from a shell, with no Python, through the magnus command.

What it has been used for

Each of these is one call with a different Hamiltonian, profile or observable.

• Beam experiments — appearance probabilities along the DUNE, T2K, Hyper-K and ESS chords, from two named sites (notebook 04).

• Atmospheric oscillograms — probability over zenith angle and energy in a single batched call (notebook 06).

• Solar neutrinos — twelve standard solar models, taken by name, and the averaged probability an experiment sees (notebook 13).

• Supernova shock fronts — where a travelling discontinuity changes the conversion probability itself (notebook 14).

• Astrophysical flavor composition, including pseudo-Dirac pairs that stay coherent after everything else has averaged (notebook 29).

• A Hamiltonian of your own — a long-range \(L_e - L_\mu\) interaction sourced by the Sun’s electrons, a cavity in the Earth’s crust, geoneutrinos, or a jet inside a collapsing star (notebook 19).

Numerical recipes gives the code for each in a few lines; Tutorial notebooks is the guided tour.

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 discretizes. 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 the same call at 1e-7	median 9e-7, worst 2e-3
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, grid and tolerances pinned	1e-12
n_jobs > 1 vs serial	exactly 0.0

The last row is the one worth reading twice: it is a bit-identity assertion rather than a tolerance, so an optimization that changed an answer would fail it rather than pass quietly. The batched scan is held to 1e-12, and only with the grid pinned, which is what isolates the batching: left to refine on its own it builds the matter profile once for the whole scan, moving the answer at the 1e-6 level. Notebook 24 measures that comparison.

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 Accuracy and diagnostics 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:

1. 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.

2. 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.

3. 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).

4. 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:

1. 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.

2. 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.

3. 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.

4. 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).

And separately from the physics, 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 answers a constant-density call exactly: the series terminates at its first term, so the evolution operator is a single exponential carrying no discretization at all, and a whole scan is one batched exponential. What NuOscProbExact offers there is not a better answer but a leaner route to it, being built for that case alone rather than carrying a general solver’s dispatch, validation and refinement ladder. The margin is narrow and runs both ways: on a 3ν constant-density scan Magνs costs 1.10 µs per energy against NuOscProbExact’s 1.44 µs batched, while at a single point it costs 33.8 µs against 19.9 µs. Against other codes carries the decision table — which of the two to reach for, case by case — together with the measured speed and accuracy behind it.

What Magνs earns its place on

The lists above are cases, not a reason. Where a closed form exists, an exact algebraic solution beats a truncated series — that is arithmetic, not a defect in anybody’s code. What is left is three axes, every one of them measured against the other codes in Against other codes:

Reach. A slab product has a floor: on a smooth profile its error bottoms out near \(2.5\times10^{-11}\) and then rises, because the round-off of composing that many matrix products outgrows what another halving of the width buys. No setting reaches below it. Magνs continues to \(2.9\times10^{-13}\).

Generality. An arbitrary \(H(t)\) — a custom Hamiltonian, a BSM term nobody has diagonalized, a profile interpolated from a simulation — needs no per-model work, because nothing in the method assumes a form for \(H\). The SU(N) closed forms stop at SU(4); Magνs has no ceiling.

Pre-packaged observables. average=True returns the phase-averaged probability a solar experiment actually measures, without resolving some 13 000 radians of phase and averaging the result yourself. Neither of the other codes offers it. And every entry point can hand back the converged evolution operator alongside the probabilities (return_evolution_operator=True), for the observables that are built from amplitudes rather than from probabilities.

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.

Performance 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 on an exponential fit or any of twelve standard solar models (Standard solar models), 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 10, with the term recursion verified term-by-term against the literature, and three integration methods. The default, Gauss-Legendre collocation integrators, reaches orders 2/4/6/8 from only 1/2/3/4 Hamiltonian evaluations per slab. Cumulative trapezoid and Simpson quadrature reach order 10 and serve 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:

• Installation & Requirements
  • Requirements
  • Installation
  • Verifying the Installation
  • File Tree
• Quick Start Guide
  • Your first probability
  • Units
  • Choosing a global fit
  • 1. Vacuum oscillations
  • 2. Matter with constant or exponential density
  • 3. The Earth (PREM) and the Sun
  • 4. Beyond the Standard Model: NSI and LIV
  • 5. Your own Hamiltonian through the Earth or the Sun
  • 6. Fully generic: any Hamiltonian, any environment

Using Magnus:

• Numerical recipes
  • One probability
  • The evolution operator, for observables built from amplitudes
  • A scan, without a loop
  • Through the Earth
  • A profile of your own
  • Phase-averaged probabilities
  • Asking for an accuracy instead of a slab count
  • Choosing a strategy, and seeing which engine answered
  • Telling it where the profile is not smooth
  • New physics
  • Writing an H_func that does not cost you a factor of five
  • Where to go next
• Tutorial notebooks
  • Start here
  • Geometry, and what experiments measure
  • New physics
  • What the method actually does
  • Where the limits are
  • Conventions worth getting right
  • Physics questions
  • Using and diagnosing the machinery
  • Against other codes, and against the parameters
  • Watching it happen
  • The paper
• Against other codes
  • The boundary, stated once
  • Which one should I use?
  • Constant density
  • PREM, three flavors
  • PREM, 3+1
  • A smooth profile: where each method runs out
  • The Sun: an observable the others do not offer
  • A supernova shock: the width of the front decides
  • Summary
  • Before comparing any two codes’ numbers
• Available Oscillation-Probability Functions
  • Vacuum
  • Matter, constant density
  • Matter, exponential density
  • Earth
  • Sun
  • Generic entry points
  • Returning the evolution operator
• Standard solar models
  • What is tabulated, and what is computed from it
  • Outside the table
  • Coherent probabilities through a table
  • Provenance
• Command-Line Calculator
  • Installation
  • Usage pattern
  • Examples
  • Choosing a propagation strategy
  • Errors are explicit rather than silent
  • Full flag reference
  • Implementation notes
• Pre-Packaged Plotting Tools
  • Why it exists
  • What it draws
  • Installation
  • A first figure
  • Adding the error subpanel
  • API conventions
  • Labels
  • Saving

How it works:

• Methodology
  • The Magnus expansion
  • Two integration methods
  • Exact unitarity from the eigendecomposition
  • Time-ordering
  • Adaptive refinement and slab placement
  • Choosing the expansion order
  • Silent vectorization and the energy-batched scan engine
  • Writing an array-capable Hamiltonian
  • Validation strategy
  • Conventions
• Magnus Expansion Terms to Any Order
  • The recursion
  • How many terms there are
  • The expansion, written out
  • Going beyond the implemented ceiling
  • Checking the implementation against the derivation
  • Choosing an order
  • API
• Adiabatic + Magnus Hybrid Strategy
  • The problem: extreme accumulated phase
  • Adiabatic transport
  • Detecting resonances without differentiating eigenvectors
  • Growing and merging non-adiabatic windows
  • Patching, and composing the pieces exactly
  • Self-certification
  • The strategy parameter
  • Validation
  • Limitations and scope
• Phase-Averaged Probabilities
  • The problem: a phase nobody can resolve
  • The averaged limit
  • The phase average
  • Coherence, and where the spread matters
  • Position-dependent Hamiltonians
  • When there is no closed form
  • Cost
  • Usage
  • Am I computing the wrong thing? strategy_info['sampling']
  • How much does the phase actually matter?
  • Limitations and scope
• Code Architecture
  • Module layout
  • The three-layer structure of magnus.oscprob
  • The layer contract: what a wrapper must not do
  • Data flow: how the Hamiltonian and potential are built
  • How to add your own wrapper
  • Where things live: a quick lookup
• Engines and dispatch
  • The engines
  • Dispatch
• Performance
  • Speed
  • How the constants were set
  • Reproducing any of this
• Accuracy and diagnostics
  • Accuracy
  • Robustness, and what each safeguard cannot do
  • Warnings: what each one means and what to do about it

Reference:

• API Reference
  • magnus.adiabatic
  • magnus.avgprob
  • magnus.cli
  • magnus.earth
  • magnus.expansionterms
  • magnus.expmkernels
  • magnus.globaldefs
  • magnus.hamiltonians
  • magnus.magnus
  • magnus.matter
  • magnus.oscprob
  • magnus.oscprobstd
  • magnus.plotting
  • magnus.solarmodels
• How to cite
  • Cite the software
  • What to say in the text
  • Citing the method
  • Related software
• References
• Changelog
  • [Unreleased]
  • [1.1.0] - 2026-09-06
  • [1.0.13] - 2026-09-06
  • [1.0.12] - 2026-09-05
  • [1.0.11] - 2026-09-05
  • [1.0.10] - 2026-09-05
  • [1.0.9] - 2026-09-05
  • [1.0.8] - 2026-09-05
  • [1.0.7] - 2026-09-05
  • [1.0.6] - 2026-09-05
  • [1.0.0] - 2026-08-11
  • [1.0.0rc1] - 2026-07-31

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.

Indices and tables

• Index

• Module Index

• Search Page