NuOscProbExact: Exact Neutrino Oscillation Probabilities

Important

Important links:

NuOscProbExact computes exact two-, three- and four-flavor neutrino oscillation probabilities for arbitrary time-independent Hamiltonians. It is a Python implementation of the method of Ohlsson and Snellman [OS00], revisited in [Bus19].

The method expands the Hamiltonian and the time-evolution operator in the bases of SU(2), SU(3) and SU(4) matrices, which yields concise, analytical, and exact closed-form expressions for the probabilities.

Four flavors brings 3+1 sterile scenarios into scope as closed four-state systems, rather than as a leak out of the three-flavor block. It is also where the method ends: Why the method stops at four explains why there is no five-flavor counterpart, and why that is a fact about algebra rather than a gap in the code.

What “exact” means here

There is no approximation beyond floating-point round-off. The evolution operator produced by the SU(2), SU(3) and SU(4) expansions is the matrix exponential \(e^{-i H_0 L}\), and the test suite checks that it is, against an independent computation:

Property	Measured agreement
\(U_3\) vs `scipy.linalg.expm`, 100 random Hermitian Hamiltonians	8e-15
\(U_2\) vs `scipy.linalg.expm`, 100 random Hermitian Hamiltonians	1.6e-15
\(U_4\) vs `scipy.linalg.expm`, 400 random Hermitian Hamiltonians	2.9e-14
Unitarity, \(U^\dagger U - \mathbb{1}\)	5e-15
Hard-coded \(d_{ijk}\) vs \(\frac{1}{4}\mathrm{Tr}(\{\lambda_i,\lambda_j\}\lambda_k)\), all 512 entries	2e-16
Exact result vs the standard vacuum oscillation formula	7e-16 (two flavors), 1e-14 (three)
Four flavors with the sterile angles off, vs `oscprob3nu`	3e-14, in vacuum and in matter alike
Three flavors vs nuSQuIDS, an independent external code	2e-15
Four flavors vs nuSQuIDS	4e-16 to 3e-10
Matter spectrum vs the Zaglauer-Schwarzer closed form	7e-16

Each figure is the worst case the corresponding test actually reaches, on the conditions that test uses — not the tolerance it asserts against, which is looser. tests/test_documented_figures.py guards the performance numbers on this page; these accuracy ones are reproduced by running the suite.

One line of context for the four-flavor row. A stiff spectrum — a 3+1 scenario with an eV-scale \(\Delta m^2_{41}\) — reaches about \(10^{-9}\) rather than \(10^{-14}\), limited by what double precision retains when the SU(4) invariants are formed rather than by the expansion itself. Both figures are far below anything an experiment can resolve, so this is a statement about the exactness claim and about error accumulating over composed slabs, not about whether the probabilities are good enough to use. Stiff spectra, and what they cost gives the mechanism, what the code does about it, and what the alternatives measured.

The last three rows are the ones an internal suite cannot supply. Everything above them compares this library against itself or against a formula transcribed from the same papers, so a convention that were wrong consistently — a mixing-matrix ordering, a sign, a unit — would pass all of it. The external checks would not, and they cover the antineutrino rule and the mass ordering explicitly. Notebook 17 works through both, including the two conventions that have to be matched first and the one residual that turns out to be ours rather than a disagreement.

When is NuOscProbExact a good fit?

Your Hamiltonian is constant, or piecewise constant. The method assumes a Hamiltonian that does not change, and in exchange gives a closed form rather than a numerical integration. A trajectory made of pieces is handled by slabs, which solves each piece exactly and multiplies the operators, and earth builds those pieces from the Preliminary Reference Earth Model.

A profile that varies smoothly over an oscillation length is the case to avoid: it can be slabbed, but the step size is then set by the oscillation rather than by the density, which for the Sun means of order \(10^4\) slabs per resonance crossing. Use a Magnus-type method there.
You want an arbitrary Hamiltonian, not a fixed scenario. The core routines take any Hermitian \(2\times2\), \(3\times3\) or \(4\times4\) matrix. Non-standard interactions, Lorentz-invariance violation, sterile states and matter effects are all just entries in that matrix; the bundled Hamiltonians in hamiltonians2nu, hamiltonians3nu and hamiltonians4nu are examples, not limitations.
You need the evolution operator, not only the probabilities. oscprob3nu.evolution_operator_3nu() returns \(U_3(L)\) itself, so it can be composed across segments or used to propagate a density matrix.
You are scanning, not evaluating one point. Every core routine takes a stack of Hamiltonians, an array of baselines, or both broadcast against each other, and returns the whole scan in one call.

What it is not

Not a solver for continuously varying Hamiltonians. See When to use Magnus instead just below, which says plainly when to reach for a different tool.
Not a five-flavor code. The expansions run to SU(4) and stop there, because the closed form does: see Why the method stops at four. Four flavors covers 3+1, which is the case people actually ask for.
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.

When to use Magnus instead

NuOscProbExact assumes the Hamiltonian is constant, or piecewise constant. Everything it is good at follows from that, and so does the one case where it is the wrong tool.

Reach for Magnus instead when the Hamiltonian varies continuously and appreciably over an oscillation length. A smoothly varying profile can always be approximated by slabs, but the step size is then set by the oscillation rather than by the density, and the slab count grows until the calculation is neither exact nor quick. Concretely:

Situation	Use this	Because
Constant density	NuOscProbExact	One closed form, no integration
Piecewise constant, tens of layers — the Earth through PREM	NuOscProbExact (`slabs`, `earth`)	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	Magnus	Slabbing needs \(\sim 10^4\) steps per resonance crossing
Genuinely open systems: decay, decoherence	Neither	Needs a Lindblad solver, not a unitary one

Notebook 14 works the solar case through and shows exactly where the wall is, rather than asserting it.

Performance

A single probability takes about 8 microseconds for three flavors and 1 for two. Scans — a curve versus baseline or energy, an oscillogram over both — are what the code mostly does, and two things make those much faster without changing any answer.

Pass arrays instead of looping. Every core routine takes a stack of Hamiltonians, an array of baselines, or both, and evaluates them in one call. That is worth roughly 20 to 90 times the equivalent Python loop, and needs no extra dependency: the expensive part of the expansion, the characteristic equation whose roots give the oscillation phases, depends on the Hamiltonian alone, so a scan over baselines solves it once instead of once per point. See Scanning: pass arrays for the three patterns.

Install Numba, if the scans are large. With pip install "nuoscprobexact[fast]" the batched paths run as compiled loops spread over the available cores; without it the NumPy path is used and the results are the same to round-off. On 2000-point scans, against the Python loop:

Scan	Loop	Arrays	Arrays + Numba
Three flavors, versus baseline	38 ms	1.8 ms	0.31 ms
Three flavors, versus energy	34 ms	1.5 ms	0.20 ms
Oscillogram, 100 x 100	197 ms	5.3 ms	0.85 ms
Two flavors, versus baseline	6.9 ms	0.07 ms	not used

Best of seven runs, interleaved, on one machine. Repeated runs vary by tens of per cent, so read these as orders of magnitude rather than constants — and if the exact figures matter to you, notebook 09 measures them on the machine that runs it.

The last row is not an omission. The backend is used only where it has been measured to win: for three flavors that is every stack size, but the two-flavor expansion reduces to a square root and a sine per element, which NumPy already does about as well as compiled code can. Below fifty thousand elements the NumPy path is kept; above it the kernel leads by about 1.3 to 1.8 times. The library chooses without being asked.

One cost runs the other way. Every entry point verifies that the Hamiltonian is Hermitian, because one that is not returns probabilities that still sum to one and so betrays nothing. Validating a stack is a pass over it — the same order of work as evaluating it, and the compiled kernel has made evaluating it fast — so it costs 1.3 to 1.8 times on a 2000-point scan and 3.2 to 5.7 times on a 200 000-point one. It is on by default anyway; where the Hamiltonians come from a construction already trusted, set oscprob3nu.CHECK_HERMITICITY to False, and likewise on the other two modules.

Methodology explains where the time goes, and what was tried and rejected.

Getting started

pip install nuoscprobexact

or, for a clone with the notebooks, the worked examples and the test suite:

git clone https://github.com/mbustama/NuOscProbExact.git
cd NuOscProbExact
pip install -e .

import numpy as np
import oscprob3nu
import hamiltonians3nu
from globaldefs import *

h_vacuum_energy_indep = \
    hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(
        S12_NO_BF, S23_NO_BF, S13_NO_BF, DCP_NO_BF, D21_NO_BF, D31_NO_BF)
h_vacuum = np.multiply(1./1.e9, h_vacuum_energy_indep)

Pee, Pem, Pet, Pme, Pmm, Pmt, Pte, Ptm, Ptt = \
    oscprob3nu.probabilities_3nu(h_vacuum, 1.3e3*CONV_KM_TO_INV_EV)

See Installation and Quickstart for the longer version.

User Guide:

Author

NuOscProbExact 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 you use NuOscProbExact in your work, please cite [Bus19]. The BibTeX entry is on the References page, and INSPIRE keeps an up-to-date record.

License

NuOscProbExact is released under the MIT License. The full text ships with the source, as LICENSE in the repository root.