Magnus

Magνs computes neutrino oscillation probabilities for two to five flavors, or for any Hermitian Hamiltonian, in vacuum, in matter, through the Earth and through the Sun. Its evolution operator is unitary to round-off at any accuracy setting, so every probability is non-negative and every row sums to one.

Installation

pip install magnuspy

Magνs requires Python 3.10 or newer. The distribution is magnuspy on PyPI, because the name magnus is taken; the import package is magnus. From a checkout, use pip install -e ., or pip install -e '.[test]' to run the tests.

Your first probability

import magnus.oscprob as oscprob
import magnus.globaldefs as gd

# 3 flavors, 1 GeV, 1300 km of vacuum
P = oscprob.osc_prob_3nu_vacuum(1.0*gd.UNIT_GEV, 1300.0*gd.UNIT_KM)
print(P[gd.NUMU][gd.NUE])        # P(nu_mu -> nu_e) = 0.0313

P[i][f] is the probability for $\nu_i \to \nu_f$, with gd.NUE, gd.NUMU, gd.NUTAU = 0, 1, 2. Energies and distances are in natural units: multiply by gd.UNIT_GEV and gd.UNIT_KM on the way in. Oscillation parameters default to the NuFIT 6.1 best fit, normal ordering.

The quick start guide takes you from here through matter, the Earth, the Sun, antineutrinos, new physics and a Hamiltonian of your own, with every example runnable as written.

More in a few lines

import numpy as np
import magnus.earth as earth

osc = gd.load_nufit_params('NuFIT 6.1')          # best-fit parameters, as sines
# Angles come back as sines (s12 = 0.556 is sin(theta12), so sin^2 = 0.309),
# phases in radians.  Leaving **osc out gives the same numbers: unset
# parameters default to this same fit.

# 200 energies through the Earth in one call: cos(theta_z) = -0.8 fixes the
# direction of the chord, and magnus.earth gives its length
E = np.linspace(0.5, 10.0, 200)*gd.UNIT_GEV
L = earth.distance_traveled_inside_earth(-0.8)*gd.UNIT_KM
P = oscprob.osc_prob_3nu_earth(E, costhz=-0.8, L=L, **osc)          # shape (200, 3, 3)

# Antineutrinos, one channel
P_bar = oscprob.osc_prob_3nu_earth(E, costhz=-0.8, L=L, nubar=True,
                                   nu_i=gd.NUMU, nu_f=gd.NUE, **osc)   # shape (200,)

# Solar neutrinos from the center of a standard solar model, phase-averaged
# over a 10% energy resolution, as a solar experiment measures them
E_sun = np.logspace(-1, np.log10(20.0), 50)*gd.UNIT_MEV
P = oscprob.osc_prob_3nu_sun(E_sun, gd.SUN_RADIUS*gd.UNIT_KM, 0.0, **osc,
                             density_profile='B16-GS98', average=True)

# The evolution operator alongside the probabilities
P, U = oscprob.osc_prob_3nu_earth(1.0*gd.UNIT_GEV, costhz=-0.8, L=L,
                                  return_evolution_operator=True, **osc)

A phase-averaged result (average=True) depends on the energy spread, average_spread, wherever some oscillation phase is small, and a PhaseAveragingWarning flags those cases. The numerical recipes page has a runnable snippet for each common task.

Or from the command line

$ magnus --flavors 3 --environment vacuum \
    --energy 1 --energy-unit GeV --baseline 1300 --baseline-unit km
Magνs 1.2.0 -- osc_prob_3nu_vacuum
E = 1 GeV, L = 1300 km

            nu_e   nu_mu  nu_tau
nu_e      0.9289  0.0085  0.0625
nu_mu     0.0313  0.3923  0.5764
nu_tau    0.0398  0.5992  0.3611

Add --scenario nsi, --flavors 5, --nu-i e --nu-f mu for one channel, or --json to pipe the result into another program. The command computes one probability at a time; scans, averages and custom Hamiltonians need Python. The CLI page is the full reference.

Why Magνs

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. An energy scan is one batched call rather than a loop; batching makes each probability one to two orders of magnitude cheaper. The cost of a probability follows how fast the density varies, not how many times the neutrino oscillates. Performance has the timings.

Accurate. Magνs propagates the evolution operator with the Magnus expansion: it exponentiates truncated integrals of the Hamiltonian over a chain of position slabs. At a tight tolerance, it agrees with an independent integration to a few parts in 10¹² at two to five flavors. Where it cannot certify its own answer, it says so.

The Against other codes page measures when to use Magνs, and when another code is cheaper.

How do I say that? It is pronounced like the name Magnus: the Greek letter ν (nu), the neutrino's symbol, stands in for the "nu" syllable.

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 measures: each interference term damped by the spread of its phase across the energy resolution, 10% by default, 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.

Examples

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, one batched energy scan per zenith angle (notebook 06).
• Solar neutrinos — standard solar models taken by name, and the averaged probability an experiment sees (notebook 13).
• Supernova shock fronts — where a traveling discontinuity changes the conversion probability itself (notebook 14).
• Astrophysical flavor composition — 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).

Every figure below is taken from the executed output of a notebook in notebooks/.

Oscillation probabilities against baseline or energy, for two to five flavors, in vacuum and in matter. notebook 03	Between two points on the Earth's surface — Fermilab to SNOLAB, Homestake, CERN and the South Pole, through PREM. notebook 04
Oscillograms across zenith angle and energy, one batched energy scan per zenith angle. notebook 06	CP violation, as biprobability ellipses traced by the CP phase. notebook 05
Five flavors: a 3+2 sterile spectrum, its fast oscillation filling the three-flavor envelope. notebook 07	A Hamiltonian of your own — here a long-range $L_e - L_\mu$ force through the Earth, against the standard curve. notebook 19
Arrangement matters: the same average density over the same path length, ordered differently, gives different probabilities. notebook 18	Phase-averaged probabilities: what remains when the oscillation is too fast for a detector to resolve. notebook 10
The averaged solar survival probability, returned directly in about 0.7 s. The green trace is the instantaneous probability that another code returns, which oscillates between 0.15 and 0.9. notebook 25	BSM against the standard curve: NSI and a sterile state on a real BS2005 solar model, with the departure below. notebook 13
The same two scenarios on a supernova shock, where the same $\varepsilon$ moves the probability by up to 0.14, against 0.014 in the Sun. notebook 14

When is Magνs not the right tool?

Magνs solves the Schrödinger equation for a Hermitian Hamiltonian fixed before the propagation. Decoherence, coupling to a bath and decay to invisible states need a non-unitary evolution; nuSQuIDS is the tool for those. Collective oscillations need a Hamiltonian that depends on the solution, and Magνs does not ship the self-consistent iteration. A feature narrower than every sampling grid goes unseen by every engine; Magνs scans the profile for such features and warns, but the scan does not catch every one. Magνs is not a flux, cross-section or detector code, a fitting framework or an event generator: it computes oscillation probabilities and stops there.

Performance

A single three-flavor probability through the Earth takes about 2 ms at the default tolerance of 10⁻³. Over 164 Earth, solar, vacuum and constant-density configurations, the median call takes 2 ms and the slowest under a second. These times are per call, on one laptop, and exclude the first call of a session. That call also loads the compiled kernels, which takes 0.1–0.3 s, or about 2 s on a machine's first run, when they compile. The timing harness is docs/dev/adversarial_batteries/timing.py.

Pass arrays instead of looping. Every wrapper takes an array of energies, of baselines or of both, and shares work across the points; the gain from batching holds at every number of flavors from two to five.

Write your H_func to accept an array of positions. It is then called once per refinement stage rather than once per quadrature node, which is several times faster, with identical output; the numerical recipes show how. A Hamiltonian that ignores its argument is detected and broadcast automatically.

n_jobs helps only where no batched engine applies: there, ten workers make a scan 2–3× faster. Where a batched engine applies, a single process is faster. See performance.

What "accurate" means here

Unitarity holds at any order and tolerance (see above). Accuracy is measured against closed forms, an independent ODE solver, expansion terms generated independently from the recursion, and populations of random profiles; the table of checks and results is in the documentation. The tolerances rtol and atol are a stopping rule, not a guarantee; see what they control.

Salient features

• Any number of flavors, any Hamiltonian: validated wrappers for 2ν, 3ν, 4ν (3+1) and 5ν (3+2), named by environment and scenario (functions), plus a generic osc_prob that takes a Hermitian Hamiltonian of any dimension.
• Vacuum, matter, Earth and Sun: constant and exponentially falling density, the Earth through PREM including chords between named sites, the Sun on an exponential fit or a tabulated standard solar model (solar models), or any profile you supply.
• The Magnus expansion to order 10, with the terms of orders 1–10 checked against terms generated independently from the recursion. The default Gauss–Legendre integrators reach orders 2, 4, 6 and 8 from 1, 2, 3 and 4 evaluations of the Hamiltonian per slab; cumulative trapezoid and Simpson quadrature reach order 10.
• Adaptive refinement to a requested tolerance, slab edges on density discontinuities, and an energy-batched engine for scans; seven engines in all, chosen from the request (engines).
• Plotting: the figures the notebooks use, one call each, from curves and small multiples to biprobability ellipses and oscillograms (plotting).

Documentation

The full documentation is at mbustama.github.io/Magnus.

Repository layout

The tree below shows the top level of the repository. The documentation also shows the package under src/. The complete listing, with a comment on every tracked file, is TREE in tests/test_file_tree.py, which generates this tree and checks it against the tracked files.

Magnus/
├── .github/                        # GitHub Actions workflows: tests, lint, notebooks, docs, publishing
├── .gitignore                      # Build, cache and generated-output artifacts
├── CHANGELOG.md                    # Version history (Keep a Changelog format)
├── CITATION.cff                    # Machine-readable citation metadata; drives GitHub's "Cite this repository"
├── LICENSE                         # GNU GPL v3 (GPL-3.0-only), the full license text
├── MANIFEST.in                     # Adds tests/conftest.py to the sdist, so that tests needing a source checkout are skipped there
├── README.md                       # Project overview; also the PyPI project description
├── docs/                           # Sphinx documentation configuration and source
├── fig/                            # Plots produced by the example notebooks
├── img/                            # Figures used by the documentation
├── notebooks/                      # Numbered Jupyter notebooks -- see docs/source/tutorials.rst
├── pyproject.toml                  # Build system, dependencies and the `magnus` console-script entry point
├── resources/                      # Paper sources and benchmark drivers, kept in the repository but not distributed
├── src/                            # The package itself -- the only thing a `pip install` delivers
├── tests/                          # Test suite (pytest; runs in CI)
└── tools/                          # Standalone utilities that are not part of the package

Continuous integration

Every push runs the test suite on Python 3.10–3.13 and lints with Ruff; the documentation is built with warnings as errors. The 29 notebooks execute too, across four parallel shards, with each shard served from cache when neither the package nor the notebook sources and data have changed. The badges at the top report those runs.

Requirements

Magνs depends on numpy (>= 1.22), scipy (>= 1.9), joblib (>= 1.2), matplotlib (>= 3.5) (so that magnus.plotting works in any installation) and numba (>= 0.59) (for the fast matrix exponential). To run the tests:

pip install -e '.[test]' && pytest tests/

Add --cov --cov-report=term-missing for coverage, configured in pyproject.toml so it measures the same thing locally and in CI.

Changelog

See CHANGELOG.md, also rendered in the docs.

How to cite

If Magνs contributed to work you are publishing, please cite it and say which version you used, since results can depend on it. The citing page has the BibTeX entry.

Mauricio Bustamante (2026). Magνs: neutrino oscillation probabilities for any Hermitian Hamiltonian, any number of flavors, and any matter profile. GitHub repository: https://github.com/mbustama/Magnus.

Methodology references:

• Sergio Blanes, Fernando Casas, José A. Oteo & José Ros (2009). The Magnus expansion and some of its applications. Physics Reports, 470(5–6), 151–238. doi:10.1016/j.physrep.2008.11.001.
• Sergio Blanes, Fernando Casas & José Ros (2000). Improved high order integrators based on the Magnus expansion. BIT Numerical Mathematics, 40(3), 434–450. doi:10.1023/A:1022311628317.

License

GNU General Public License v3.0 only (GPL-3.0-only); the full text is in LICENSE. You may use, study, modify and redistribute it, including commercially, provided derivative works carry the same license with source available.

Author

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