python

Magnus: neutrino oscillation probabilities for any Hamiltonian, in any number of flavors

Standard

Documentation: https://mbustama.github.io/Magnus/

Notebooks: https://mbustama.github.io/Magnus/tutorials.html

Source code: https://github.com/mbustama/Magnus

PyPI: https://pypi.org/project/magnuspy/

Oscillations in matter are usually computed by pretending the matter holds still: chop the path into slabs of constant density, solve each one exactly, multiply. It is a good trick, and for a neutrino crossing the Earth it is nearly the truth. For the Sun, for a supernova, for anything with a front moving through it, the pretending is where the physics is.

Magνs is my new package for that case. It computes oscillation probabilities for an arbitrary number of flavors, under any Hamiltonian, whether or not it depends on time. Rather than stepping the Schrödinger equation forward, it integrates the Hamiltonian across each slab and exponentiates the result — the Magnus expansion. The useful consequence is structural: however hard the expansion is truncated, the truncation stays inside the Lie algebra, so the evolution operator is exactly unitary by construction. Probabilities come out non-negative and summing to one at machine precision, at every accuracy setting, not only the expensive ones.

What it makes easy

  • A density that actually changes — a tabulated solar model, a supernova profile from a simulation, a shock front.
  • An accuracy instead of a slab count — ask for a tolerance and the refinement finds the resolution for you.
  • Any number of flavors — two through five, sterile states (3+1, 3+2) included.
  • The phase-averaged probability, in closed form — the quantity a solar experiment actually measures, returned directly rather than reconstructed by averaging a scan.
  • Neutrinos through the Earth — you give a direction, not a profile.

Install with:

pip install magnuspy

The import package is magnus. Python 3.10 or newer, GPL-3.0-only, and 27 example notebooks that are executed in CI, so the numbers printed in them cannot quietly drift.

A minimal example

A 10 GeV muon neutrino arriving from below, crossing the Earth:

import magnus.oscprob as oscprob
import magnus.globaldefs as gd
import magnus.earth as earth
# A 10 GeV muon neutrino crossing the Earth, arriving from below.
costhz = -0.4
L = earth.distance_traveled_inside_earth(costhz)*gd.UNIT_KM   # chord: km -> eV^-1
P = oscprob.osc_prob_3nu_earth(10.0*gd.UNIT_GEV, costhz=costhz, L=L)
print('Pme = %.5f, Pmm = %.5f, Pmt = %.5f'
      % (P[gd.NUMU][gd.NUE], P[gd.NUMU][gd.NUMU], P[gd.NUMU][gd.NUTAU]))

which prints

Pme = 0.10273, Pmm = 0.01165, Pmt = 0.88562

The oscillation parameters are never mentioned, because they default to the current global fit (NuFIT 6.1). Neither is the Earth’s density: the wrapper walks the PREM profile for the chord that the zenith angle picks out. What you supply is a direction and an energy.

How it differs from NuOscProbExact

Both packages are mine, and they are not competitors — the honest way to put it is as a boundary rather than a winner. NuOscProbExact solves each constant-density slab in closed form, through SU(2), SU(3) and SU(4) expansions. It is exact up to round-off, there is no slab count to argue about, and it is fast. Where a closed form exists and the accumulated phase is large, use the closed form — including the Earth through PREM at three flavors, where NuOscProbExact is about 20× cheaper per call than Magνs. That is its home ground.

Magνs integrates across each slab instead, which is exactly what lets it follow a density that changes while the neutrino is still inside it. Reach for it when the profile varies continuously, when you would rather ask for an accuracy than guess a resolution, when you need five flavors — the SU(N) closed forms stop at SU(4) — or when the Hamiltonian is something nobody has diagonalized.

Two measurements make the difference concrete. On a smooth profile, multiplying constant-density slabs together bottoms out near 2.5 × 10⁻¹¹ and then gets worse: past roughly 16 000 slabs, round-off from composing that many matrix products outgrows anything a finer grid buys. Magνs keeps going to 2.9 × 10⁻¹³. And for the Sun, it returns 40 phase-averaged energies in about 0.7 s, where nuSQuIDS needs on the order of ten minutes merely to reach the solver tolerance at which its output is a probability at all.

The two packages deliberately share conventions, units and parameter defaults, so a calculation can be moved from one to the other and used as a cross-check on itself. I do this often.

Some physics it’s meant for

Each of these has a worked notebook in the tutorial gallery, and there is a systematic comparison against other codes in the documentation.

  • Long-baseline beams and atmospheric neutrinos through the Earth.
  • Solar neutrinos on a real tabulated model, including the averaged observable.
  • Supernova neutrinos across a shock front — where the width of the front is what decides which method the problem belongs to.
  • Sterile states (3+1 and 3+2), non-standard interactions, Lorentz-invariance violation.
  • A Hamiltonian of your own — the escape hatch, and the reason the rest of the list is not exhaustive.

NuOscProbExact 1.11.0: exact neutrino oscillation probabilities for two, three, and now four flavors — on PyPI!

Standard

Documentation: https://mbustama.github.io/NuOscProbExact/

Notebooks: https://github.com/mbustama/NuOscProbExact/tree/main/notebooks

Source code: https://github.com/mbustama/NuOscProbExact

PyPI: https://pypi.org/project/nuoscprobexact/

Neutrino oscillation probabilities look like a solved problem until the Hamiltonian stops being the textbook one. There is a closed form in vacuum, and another in two-flavor constant-density matter, and after that it runs out. Add a realistic density profile, a non-standard interaction, a Lorentz-violating operator or a fourth state, and the usual options are to expand in a parameter that may not be small, to diagonalise the Hamiltonian numerically at every point of a grid, or to integrate the Schrödinger equation and hope the step was short enough.

NuOscProbExact does none of these. The Hamiltonian and the time-evolution operator are expanded in the basis of SU(2), SU(3) or SU(4) matrices, and exponentiated analytically using the closure of the algebra. What comes out is a closed form for the probabilities in terms of two or three invariants of the Hamiltonian: no eigenvectors, no numerical diagonalisation, no integration, and no approximation beyond floating-point round-off. Any Hermitian 2×2, 3×3 or 4×4 matrix will do, whatever physics put it there. The method is Ohlsson and Snellman’s, revisited in the 2019 paper; the extension to SU(4) is new since then.

Four flavors is where it stops, for a reason external to neutrino physics: at five, the eigenvalues stop being expressible in radicals — Abel–Ruffini — and the closed form ends with them. Four is exactly what 3+1 needs.

If your reference point is a numerical evolution code, the output is the same object — probabilities, or the evolution operator itself — but with no step size to converge and no dependency beyond numpy. I wrote the first version because adding a model should mean writing down a matrix, not writing a solver. Seven years on it does rather more than it did, and version 1.11.0 now installs with pip.

What it makes easy

  • Any model, without a new code path. Vacuum, matter, non-standard interactions, Lorentz-invariance violation, sterile states: none of these is a special case in the code. Each is a different matrix handed to the same routine.
  • A whole scan in one call. A curve against energy or baseline, or a full oscillogram over both, with no Python loop.
  • Neutrinos through the Earth. Along a zenith angle, or between two of fifteen named sites — CERN–Gran Sasso, Fermilab–Homestake, Tokai–Kamioka and others — with the density profile taken from PREM. Or through any hand-built sequence of slabs of arbitrary width and density, each solved exactly and the operators multiplied.
  • 3+1 sterile scenarios at four flavors, as closed four-state systems rather than as leakage out of the three-flavor block.
  • The evolution operator itself, not only probabilities, for composing across segments or propagating a density matrix.

Install with:

pip install nuoscprobexact

and, for the optional compiled backend:

pip install "nuoscprobexact[fast]"

MIT licence. Python 3.9 or newer, tested on 3.9 through 3.13, and numpy is the only required dependency: no compiler, nothing to build.

Main features

Accuracy

  • Cross-checked against independently written code. The probabilities agree with nuSQuIDS to round-off once conventions are matched, and the matter spectrum agrees with the Zaglauer–Schwarzer closed form.
  • A regression suite of 596 tests, at 100% coverage.

Speed

  • Arrays instead of loops. Stacks of Hamiltonians, arrays of baselines, or both broadcast together: roughly 20x to 90x faster than the equivalent Python loop, with no extra dependency.
  • An optional compiled backend. With numba installed, the batched paths run as compiled kernels — answers identical to round-off, used only where it has been measured to win.

What you need to provide

You provide a Hermitian matrix and a baseline, built by hand or with one of the bundled builders for the standard cases. The package returns the exact evolution operator and the probabilities that follow from it, and stops there: no fluxes, no cross sections, no detector response, no fitting. One assumption is worth stating plainly: the Hamiltonian is taken to be time-independent over each segment. A smoothly varying profile can be slabbed, but the step size is then set by the oscillation rather than by the density, and there the documentation points at a Magnus-type method instead.

A minimal example

import numpy as np
import oscprob3nu
import hamiltonians3nu
from globaldefs import *
# The vacuum Hamiltonian, without the 1/E factor, so an energy scan
# computes it once
h_vac = hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(
    S12_NO_BF, S23_NO_BF, S13_NO_BF, DCP_NO_BF, D21_NO_BF, D31_NO_BF)
# One probability, at 1 GeV over 1300 km
Pee, Pem, Pet, Pme, Pmm, Pmt, Pte, Ptm, Ptt = oscprob3nu.probabilities_3nu(
    np.asarray(h_vac)/1.e9, 1300.0*CONV_KM_TO_INV_EV)
print('Pee = %.5f, Pem = %.5f, Pet = %.5f' % (Pee, Pem, Pet))

prints

Pee = 0.92768, Pem = 0.01432, Pet = 0.05800


(The quickstart goes through this more slowly; the recipes page has short, complete snippets for most of what the package can compute, each with real output.)

What is new since 1.0.0

The original announcement, in May 2019, covered the method paper and the code on GitHub. The method has not changed since; almost everything around it has. Whole-array scans arrived in 1.2.0, the optional numba backend in 1.6.0, the Earth and arbitrary layered matter in 1.8.0, and four flavors — oscprob4nu, hamiltonians4nu, the SU(4) expansion — in 1.9.0. In 1.11.0, a non-Hermitian Hamiltonian raises instead of quietly returning probabilities that still sum to one. Along the way, a documentation site and eighteen worked notebooks.

One thing to flag for anyone still running 1.0.0: two fixes in 1.1.0 changed numbers, both at two flavors. probabilities_2nu was missing the h2 contribution to P(νe → νμ), which is invisible whenever the off-diagonal entry is real — so vacuum and constant-density matter were unaffected, but a CP-violating two-flavor Hamiltonian was not — and the two-flavor vacuum Hamiltonian used the opposite sign convention, invisible in vacuum and consequential once a matter potential is added. Three-flavor results from 1.0.0 are unaffected; the changelog has the details.

Some physics it’s meant for

Long-baseline accelerator experiments. DUNE-like baselines through the crust, with matter effects and the CP phase, and the bi-probability ellipses above, which show how matter pushes neutrinos and antineutrinos off the CP-symmetric diagonal. Notebook 05.

Atmospheric neutrinos through the Earth. Zenith-angle scans with the density taken from PREM rather than from a single average value, which is where the mass ordering and the θ23 octant separate. Notebooks 07 and 12.

Sterile neutrinos, 3+1. Short-baseline appearance and disappearance, and the sterile matter resonance through the Earth, with the four states treated as a closed, unitary system rather than as leakage out of the three-flavor block. Notebook 16.

Non-standard interactions and Lorentz-invariance violation. Both enter as extra Hermitian terms added to the vacuum Hamiltonian, so a new operator means a new matrix, not a new solver. Notebook 03.

Unusual matter profiles. Castle-wall and other hand-built profiles, where the arrangement of the matter changes the answer even when the mean density does not — a parametric-resonance effect no constant-density approximation can show. Notebook 08.

PyFC: a Python framework for Feldman-Cousins confidence intervals

Standard

Links

Setting a confidence interval sounds like a solved problem, until you actually have to do it. Deciding whether to quote a two-sided interval or a one-sided upper limit after looking at the data — the classic “flip-flopping” — breaks coverage, and near a physical boundary (a mass, a cross section, a flavor fraction) the usual construction can return an empty or unphysical interval. Feldman and Cousins’ 1998 unified approach fixed this by using the likelihood-ratio ordering principle to let the data decide the interval’s shape automatically, with guaranteed coverage throughout. It has become the standard approach for small-signal counting experiments across particle and astroparticle physics — but implementing it correctly, especially in more than one dimension and with realistic, finite-statistics systematics, is a substantial piece of software to get right.

That’s what I built PyFC to handle — it supports binned and unbinned likelihoods natively, in the same framework, rather than treating one as a special case of the other. Given a likelihood, it runs the Monte Carlo pseudo-experiments needed to empirically calibrate the profile-likelihood-ratio distribution, and returns exact-coverage Feldman-Cousins intervals and contours — in one or several dimensions, on a laptop or on a cluster. I kept needing exactly this — correct-coverage contours for neutrino mixing-parameter extractions and HNL coupling limits in my own papers — and rewriting the toy-generation machinery from scratch for each project wasn’t sustainable.

If your reference point is ROOT/RooStats’ FeldmanCousins calculator: PyFC does the same statistical job with no ROOT dependency — plain NumPy/SciPy arrays and a pip install, so it drops straight into a normal Python analysis.

Main features

PyFC is built for speed, flexibility, and accuracy — none traded off for the others.

Speed

  • Parallel processing by default. Binned fits use Numba/NumPy C-extensions with GIL-bypassing thread pooling; unbinned fits parallelize across CPU cores via Python’s multiprocessing (ProcessPoolExecutor). Toy generation — normally the bottleneck in any Feldman-Cousins analysis — scales with however many cores you give it.
  • No wasted compute. Adaptive early-stopping (adaptive_toys) and batched dispatch (toy_batch_size) stop generating toys once the result is statistically clear, instead of running a fixed, oversized number every time.
  • Resumable. Checkpointing saves state after each parameter/pair scan, so a run killed by a cluster walltime limit picks up where it left off instead of restarting from zero.

Flexibility

  • Native binned and unbinned support, in one framework. Poisson binned data — including arbitrary N-dimensional, even non-contiguous, histograms — and Extended Unbinned Maximum Likelihood models are both first-class, not one bolted onto the other.
  • A choice of optimizer. Gradient-based SciPy (L-BFGS-B), nested sampling via UltraNest, brute-force grid scans, or a Hybrid strategy that uses UltraNest’s robustness for the one global fit and SciPy’s speed for the many conditional fits that follow.
  • Physical constraints between parameters, not just per-parameter boxes. bounds_func and constraints express joint or simplex conditions (e.g. a + b ≤ 1) that an independent [lo, hi] range per parameter can’t capture.

Accuracy

  • Coverage that’s actually exact. The test-statistic distribution is derived empirically from generated toys rather than assumed from Wilks’ theorem, and disconnected accepted regions are reported as separate intervals instead of being silently merged into one.
  • Finite-Monte-Carlo corrections for limited simulation statistics, so a small MC background template doesn’t bias the result.

Also: optional 2D grid sparsification for faster contours, and an interactive CLI for generating run configurations.

Open source (GPLv3), Python 3.8+.

What you need to provide

You bring the physics, PyFC handles the statistics:

  • A physics mapper — for binned data, a compute_rates_func that maps your free parameters to expected bin counts and variances; for unbinned data, that plus pdf_components (one function per signal/background component) and a generate_toy_func for parametric bootstrapping.
  • Your data — a histogram (any shape, any dimension) or an array of events.
  • grids — the scan grid for each free parameter, plus bounds_func/constraints for any joint constraints between them.
  • A config — confidence level, toy settings, checkpointing, and the optimizer strategy (grid, scipy, ultranest, or hybrid), built interactively with the CLI or as a JSON file.

Call compute_fc_intervals and PyFC takes care of the rest: toy generation, parallelization, checkpointing, and exact-coverage intervals out the other end — no hand-rolled Monte Carlo loop required.

A minimal example

For a binned Poisson analysis with a signal and a background component:

import numpy as np
from numba import njit
from pyfc import compute_fc_intervals
# A fixed expected-shape array, captured via closure (not itself a fit parameter)
expected_shape = np.array([0.1, 0.5, 2.0, 5.0])
@njit(fastmath=True, nogil=True)
def compute_rates(params):
    """Map parameters -> expected bin counts (mu) and variances (sigma2)."""
    mu = params[0] * params[1] * expected_shape + params[2]   # signal (rate x efficiency) + background
    sigma2 = mu                                                # Poisson variance
    return mu, sigma2
observed_counts = np.array([20, 7, 2, 0])   # data — any shape, not just 1D
grids = [
    np.linspace(1e-9, 1e-7, 20),   # param_1
    np.linspace(2.0, 3.0, 15),     # param_2
    np.linspace(0.8, 1.2, 10),     # param_3
]
# results: exact-coverage FC intervals/contours for param_1-3, plus the
# underlying profile-likelihood scan and toy statistics behind them
results = compute_fc_intervals(
    data=observed_counts,
    grids=grids,
    compute_rates_func=compute_rates,
    **config,   # confidence level, toy settings, strategy — see Configuration
)

(See the Quick Start Guide for the full binned and unbinned walkthroughs, with real output.)

Install with:

pip install PyFeldmanCousins

Some physics it’s meant for

Feldman-Cousins intervals matter most exactly where physics tends to live: near a boundary, with a marginal signal, where Gaussian/Wilks assumptions don’t hold. A few non-trivial cases PyFC is a natural fit for:

  • Astrophysical neutrino flavor composition. The flavor fractions (f_e, f_μ, f_τ) of the diffuse high-energy neutrino flux live on a simplex, f_e + f_μ + f_τ = 1 — precisely the joint-constraint case that bounds_func/constraints were built for. Feldman-Cousins contours on the flavor triangle are the standard way to compare an IceCube-like measurement against competing production scenarios (pion decay, muon-damped decay, neutron decay) and BSM flavor physics. Worked through hands-on in tutorial notebook 06.
  • Heavy Neutral Lepton searches at beam-dump experiments. A null result from a fixed-target experiment like SHiP constrains the active-sterile mixing |U|² as a function of HNL mass — an upper limit sitting right against the |U|² ≥ 0 boundary, often with disconnected sensitivity islands where the excluded region reappears at higher mass — see tutorial notebook 04 for how PyFC reports these correctly.
  • Dark matter direct-detection cross-section limits. The textbook case: a handful of candidate events over an uncertain background, where the interval needs to turn two-sided smoothly if a real excess ever shows up, instead of flip-flopping into an artificially tight one. The quickstart tutorial notebook works through essentially this case end to end.
  • Accelerator and reactor neutrino oscillations. Joint 2D regions on (Δm², sin²2θ) — this is literally the original 1998 Feldman-Cousins worked example, and still a good stress test for comparing the grid, SciPy, and Hybrid strategies — see tutorial notebook 05.
  • Neutrinoless double-beta decay. Binned searches for a peak near the Q-value, where the finite-MC correction matters because the simulated background templates are themselves statistics-limited — see tutorial notebook 04 for the correction on vs. off.