Code to compute exact two-, three- and four-neutrino oscillation probabilities using SU(2), SU(3) and SU(4) expansions
Note: The oscillation probabilities are computed exactly, with no approximation beyond floating-point round-off. A regression test suite lives in
tests/and can be run withpytest, and the results are cross-checked against nuSQuIDS, an independent external code — see notebook 17.
Every figure below is produced by a notebook in notebooks/, and the link under each one goes to the code that drew it. The documentation collects the same material, with runnable snippets, on its numerical recipes page.
 Oscillation probabilities against baseline or energy, for two, three or four flavors. notebook 02 |
 Matter, NSI and Lorentz-invariance violation — each just a different Hermitian matrix. notebook 03 |
 Oscillograms over energy and baseline: 57 600 probabilities in a single call. notebook 04 |
 CP violation, as bi-probability ellipses in vacuum and in matter. notebook 05 |
 The Earth's density, from the Preliminary Reference Earth Model. notebook 06 |
 Neutrinos through the Earth, in energy and zenith angle, or between two named sites. notebook 07 |
 Arbitrary matter profiles — castle walls and worse, exactly. notebook 08 |
 Mass ordering and the θ₂₃ octant, separated by matter through the Earth. notebook 12 |
 Four flavors: a 3+1 sterile state, resolved at a short baseline. notebook 16 |
 The sterile matter resonance through the Earth, in energy and zenith angle. notebook 16 |
 A Hamiltonian of your own, carried through a varying profile — here a long-range force through the Earth. notebook 20 |
NuOscProbExact is a Python implementation of the method developed by Ohlsson & Snellman to compute exact neutrino oscillation probabilities for arbitrary time-independent Hamiltonians. The method was revisited and the code presented in the paper NuOscProbExact: a general-purpose code to compute exact two-flavor and three-flavor neutrino oscillation probabilities (arXiv:1904.12391), by Mauricio Bustamante.
The paper covers two and three flavors; the code has since been extended to four, through the SU(4) algebra, which brings 3+1 sterile scenarios into scope. Four is where the closed form ends — see why.
The method relies on expansions of the Hamiltonian and time-evolution operators in terms of SU(2), SU(3) and SU(4) matrices in order to obtain concise, analytical, and exact expressions for the probabilities, that are also easy to implement and evaluate. For details of the method, see the paper above; the four-flavor extension is documented in the methodology page.
- Exact probabilities for any Hermitian 2×2, 3×3 or 4×4 Hamiltonian. There is no approximation beyond floating-point round-off, and the probabilities agree with the independently written nuSQuIDS to round-off once conventions are matched (notebook 17). Oscillations in vacuum, in matter, with non-standard interactions, in a Lorentz invariance-violating background and with sterile states are not special cases in the code — each is a different matrix handed to the same routine.
- Four flavors, for 3+1 sterile scenarios.
oscprob4nucarries the same closed-form treatment to SU(4), which is the last place it reaches: at five flavors the eigenvalues stop being expressible in radicals, and that is a theorem rather than a missing feature. A 3+1 system is closed and unitary over all four states, so it sits squarely inside the method's assumptions rather than "leaking" out of a three-flavor block. - The evolution operator itself, not only the probabilities, so it can be composed across segments or used to propagate a density matrix.
- Whole scans in one call. Every core routine accepts a stack of Hamiltonians, an array of baselines, or both broadcast against each other, which is tens of times faster than the equivalent Python loop and gives identical results.
- Piecewise-constant matter.
slabspropagates across a sequence of adjacent slabs of arbitrary width and density, solving each exactly and multiplying the operators, at two, three or four flavors. - The Earth.
earthbuilds those slabs from the Preliminary Reference Earth Model, and computes probabilities along a given zenith angle or between two of fifteen predefined locations. A 3+1 crossing is included: the sterile state does not feel the neutral-current potential, so that potential stops cancelling andearthbuilds it per slab. PREM fixes the density but not the composition, so the electron fraction is the caller's: it defaults to one half throughout, andearth.electron_fraction_premsupplies the value of each layer instead — core, mantle, crust and ocean. That default will change to the layered values in a future release. - A compiled backend, installed by default. Large batched calls run on compiled kernels, and so does the whole slab composition behind
slabsandearth— an Earth crossing is 7–14× quicker depending on the flavor count.numbacomes with the package as of 1.13.0; where it is unavailable or switched off, the NumPy path is used and the answers are the same to round-off.
- Hamiltonians that vary continuously along the trajectory. See When to use Magnus instead below — this is the one case where a different tool is the right answer, and it is worth knowing before you start.
- More than four flavors. The expansions run to SU(4) and stop, because the closed form does: solving for the eigenvalues means solving the characteristic polynomial in radicals, and at degree five Abel–Ruffini says that cannot be done. Four flavors covers 3+1, which is the case people actually ask for.
- Neutrino production, cross sections, fluxes or detector response. This computes oscillation probabilities and nothing downstream of them.
- Fitting or statistics. There is no likelihood machinery here; the probabilities are meant to be handed to whatever does that.
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.
Use Magnus instead when the Hamiltonian varies continuously and appreciably over an oscillation length. A smoothly varying profile can always be approximated by slabs, but then the step size is set by the oscillation rather than by the density, and the slab count grows until the calculation is neither exact nor quick.
| Situation | Use | 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 ~104 steps per resonance crossing |
| 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.
NuOscProbExact was developed by Mauricio Bustamante. If you use it in your work, please follow the directions on Citing.
NuOscProbExact is fully written in Python 3. It uses standard modules that are available, sometimes by default, as part of most Python installations, either stand-alone or via Anaconda. numpy and numba are installed with the package, so the first four rows need nothing at all; where a row names an extra, that is the one to install. The commands are under Installation below, and are not repeated here.
| To do this | You need | Extra |
|---|---|---|
Compute probabilities (oscprob2nu.py, oscprob3nu.py, oscprob4nu.py) |
numpy, cmath |
— |
Use the bundled sample Hamiltonians (hamiltonians2nu.py, hamiltonians3nu.py, hamiltonians4nu.py) |
numpy, cmath, copy |
— |
Propagate through layered matter or the Earth (slabs.py, earth.py) |
numpy |
— |
| Go faster on large scans | numba |
— (installed by default) |
Run the notebooks (notebooks/) |
matplotlib, Jupyter |
notebooks |
Run the regression suite (tests/) |
pytest, scipy, coverage, pytest-cov |
test |
| Build the documentation | Sphinx and friends | docs |
numpy and numba are installed for you; nothing else is needed to compute a probability. scipy is used by the test suite alone, to cross-check the evolution operator against an independent matrix exponential; the library itself never imports it.
numba became a real dependency in 1.13.0 rather than an optional extra, because it is worth roughly 1.5x to 20x on large scans and there is no reason for that to depend on reading far enough down this page. The NumPy path remains, is still tested on every push, and is what runs if numba is unavailable or switched off — the results are identical to round-off either way. One consequence to know about: numba declares its own ceiling on numpy (0.66.0 requires numpy<2.5), so installing into an environment that holds a newer numpy than that will downgrade it. If you need the newest numpy more than you need the speed, install with --no-deps and add numpy yourself, or set fastkernels.USE_NUMBA = False and ignore the compiled path.
NuOscProbExact is pure Python: there is nothing to compile or link.
Python version: The code requires Python 3.9 or newer, and every release is tested on 3.9, 3.10, 3.11, 3.12, and 3.13. The floor comes from
numpy.broadcast_shapes, which the batched paths use and which arrived in NumPy 1.20; 3.9 is also the oldest version for which the optionalnumbabackend still has a wheel.
pip install nuoscprobexactThat is the whole installation, and it includes the compiled backend: numpy and numba both arrive, and large scans run on compiled kernels without anything further being asked for.
The optional extras add what each task needs, and can be combined:
pip install "nuoscprobexact[fast]" # numba; now a base dependency, so this is a no-op
pip install "nuoscprobexact[notebooks]" # Jupyter, matplotlib and scipy, for notebooks/
pip install "nuoscprobexact[test]" # pytest, scipy and coverage, to run the suite
pip install "nuoscprobexact[docs]" # Sphinx and friends, to build the documentationThen, in your own code:
import numpy as np
import oscprob3nu
import hamiltonians3nu
import globaldefs as gd
h_vacuum = hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(
gd.S12_NO_BF, gd.S23_NO_BF, gd.S13_NO_BF, gd.DCP_NO_BF,
gd.D21_NO_BF, gd.D31_NO_BF)
prob = oscprob3nu.probabilities_3nu(
np.asarray(h_vacuum)/1.e9, 1300.0*gd.CONV_KM_TO_INV_EV)The modules are installed under their bare names --- oscprob2nu, oscprob3nu, oscprob4nu, hamiltonians2nu, hamiltonians3nu, hamiltonians4nu, globaldefs, fastkernels, slabs, earth --- which is the same way the paper and the worked examples refer to them.
Install from a clone if you want the notebooks, the worked examples from the paper, the regression suite, or a version that is not yet released:
git clone https://github.com/mbustama/NuOscProbExact.git
cd NuOscProbExact
pip install -e .-e installs in editable mode, so edits to src/ take effect without reinstalling. The extras work the same way, for example pip install -e ".[fast,test]".
A clone gives you the following file structure:
NuOscProbExact/
├── .github/ # Continuous integration (GitHub Actions)
│ └── workflows/
│ ├── tests.yml # The suite: five Pythons, all three backends
│ ├── lint.yml # ruff, and the docs build under -W
│ ├── pages.yml # Builds and deploys the docs to GitHub Pages
│ └── publish.yml # Publishes to PyPI on a GitHub Release
├── .gitignore # Build, cache, and generated-output artefacts
├── CHANGELOG.md # Notable changes, rendered as a docs page
├── RELEASE_NOTES.md # What changed since the published version
├── LICENSE # MIT license
├── README.md # Project overview and worked examples
├── pyproject.toml # Packaging metadata and pytest configuration
├── MANIFEST.in # Keeps the test suite out of the source distribution
├── examples/ # Runnable scripts, one per scenario, linked from README.md
│ ├── example_2nu_trivial.py # Two-flavor, arbitrary Hamiltonian
│ ├── example_2nu_vacuum.py # Two-flavor, oscillations in vacuum
│ ├── example_2nu_vacuum_coeffs.py # Two-flavor, expansion coefficients
│ ├── example_3nu_trivial.py # Three-flavor, arbitrary Hamiltonian
│ ├── example_3nu_vacuum.py # Three-flavor, oscillations in vacuum
│ ├── example_3nu_vacuum_coeffs.py # Three-flavor, expansion coefficients
│ ├── example_3nu_matter.py # Three-flavor, oscillations in matter
│ ├── example_3nu_nsi.py # Three-flavor, matter with NSI
│ └── example_3nu_liv.py # Three-flavor, LIV background
├── docs/ # Sphinx documentation
│ ├── Makefile # `make html` on Linux and macOS
│ ├── make.bat # `make html` on Windows
│ ├── requirements.txt # Documentation-only dependencies
│ └── source/
│ ├── conf.py # Sphinx configuration
│ ├── index.rst # Landing page
│ ├── installation.rst # Requirements, installation, file tree
│ ├── quickstart.rst # Shortest path to a probability
│ ├── recipes.rst # Numerical recipes, with pre-generated figures
│ ├── notebooks.rst # The twenty tutorial notebooks, in reading order
│ ├── methodology.rst # The SU(2), SU(3) and SU(4) expansions
│ ├── functions.rst # API reference, from the docstrings
│ ├── references.rst # Bibliography
│ ├── refs.bib # BibTeX entries for the bibliography
│ ├── changelog.rst # Includes the root CHANGELOG.md
│ └── _static/
│ ├── nuoscprobexact_logo.png
│ └── slabs_composition.svg # How slabs compose, drawn for quickstart.rst
├── img/ # Figures from earlier versions of README.md
│ ├── anim_cp.gif # The demonstration reel, one scene per file
│ ├── anim_earth.gif
│ ├── anim_slabs.gif
│ ├── anim_sterile.gif
│ ├── prob_3nu_vacuum_vs_baseline_ee_em_et.png
│ ├── prob_3nu_vacuum_vs_energy_ee_em_et.png
│ └── gallery/ # Figures lifted from the notebooks, shown in README.md
│ ├── gallery_anim_cp.png
│ ├── gallery_anim_earth.png
│ ├── gallery_anim_slabs.png
│ ├── gallery_anim_sterile.png
│ ├── gallery_biprobability.png
│ ├── gallery_earth.png
│ ├── gallery_long_range.png
│ ├── gallery_matter.png
│ ├── gallery_ordering.png
│ ├── gallery_oscillogram.png
│ ├── gallery_prem.png
│ ├── gallery_profiles.png
│ ├── gallery_sterile.png
│ ├── gallery_sterile_earth.png
│ └── gallery_vacuum.png
├── notebooks/ # Worked examples, with their figures stored inline
│ ├── 01_basics.ipynb # Units, one probability, and broadcasting
│ ├── 02_vacuum_oscillations.ipynb # Against baseline and against energy
│ ├── 03_matter_nsi_liv.ipynb # Constant-density matter, NSI, and LIV
│ ├── 04_oscillogram.ipynb # Energy-baseline maps in one call
│ ├── 05_biprobability.ipynb # CP ellipses, in vacuum and in matter
│ ├── 06_earth_and_prem.ipynb # PREM, chord geometry, and slabs
│ ├── 07_earth_probabilities.ipynb # Through the Earth, and between sites
│ ├── 08_unusual_density_profiles.ipynb # Castle-wall and other hand-built profiles
│ ├── 09_performance.ipynb # Looping vs broadcasting, and the backend
│ ├── 10_paper_figures.ipynb # The two figures from arXiv:1904.12391
│ ├── 11_exact_vs_approximations.ipynb # Where the textbook formulas break down
│ ├── 12_ordering_and_octant.ipynb # Normal vs inverted, and the 23 octant
│ ├── 13_antineutrinos.ipynb # Conjugate and flip, and two ways to slip
│ ├── 14_solar_and_adiabatic_msw.ipynb # The MSW resonance, and the cost wall
│ ├── 15_numerical_edge_cases.ipynb # Degeneracies, and what does not go NaN
│ ├── 16_four_neutrinos.ipynb # A 3+1 sterile state, through the SU(4) expansion
│ ├── 17_cross_checks.ipynb # Corroboration from nuSQuIDS and Zaglauer-Schwarzer
│ ├── 18_evolution_operator.ipynb # The operator, and the SU(n) coefficients
│ ├── 19_animations.ipynb # Four animated scenes, as stills; the reel they came from
│ ├── 20_arbitrary_hamiltonian.ipynb # A long-range force, carried through three profiles
│ ├── make_notebooks.py # Generates and executes all of the above
│ └── coastlines_crude.json # Coastlines for the Earth cutaway, from GSHHS
├── src/ # The library
│ ├── oscprob2nu.py # Two-flavor probabilities, SU(2) expansion
│ ├── oscprob3nu.py # Three-flavor probabilities, SU(3) expansion
│ ├── oscprob4nu.py # Four-flavor probabilities, SU(4) expansion
│ ├── hamiltonians2nu.py # Example two-flavor Hamiltonians
│ ├── hamiltonians3nu.py # Example three-flavor Hamiltonians
│ ├── hamiltonians4nu.py # Example four-flavor (3+1) Hamiltonians
│ ├── globaldefs.py # Physical constants and unit conversions
│ ├── fastkernels.py # Numba kernels, with a NumPy fallback
│ ├── slabs.py # Propagation across adjacent slabs
│ └── earth.py # PREM, chord geometry, and Earth crossings
├── tools/ # Scripts that are not part of the package
│ └── make_demo_video.py # Joins and shrinks the clips notebook 19 renders
├── resources/ # Things the repository carries but does not install
│ ├── paper/ # Source of the paper that documents this library
│ │ ├── README.md # How to build it, and what each file is
│ │ ├── HANDOVER-magnus-fig11.md # What another paper needs to redraw Figure 11 with its own code on it
│ │ ├── main.tex # The paper; BibTeX, no inlined bibliography
│ │ ├── make_versions.py # Derives the clean and the diffed versions
│ │ ├── baseline_cpc_v1.tex # The published version, for the diff
│ │ ├── refs.bib # The bibliography, read by BibTeX at build
│ │ ├── elsarticle.cls
│ │ ├── elsarticle-num.bst
│ │ └── figs/ # Every figure, fourteen of them from notebook 10
│ │ ├── validation.pdf
│ │ ├── prob_2nu_vs_energy_compare.pdf
│ │ ├── prob_3nu_vs_energy_compare.pdf
│ │ ├── slabs_composition.pdf
│ │ ├── earth_oscillogram.pdf
│ │ ├── density_arrangement.pdf
│ │ ├── sterile_earth_oscillogram.pdf
│ │ ├── architecture.pdf
│ │ ├── performance.pdf
│ │ ├── exact_vs_approximations.pdf
│ │ ├── speed_accuracy.pdf
│ │ ├── scan_sensitivity.pdf
│ │ ├── prem_speed_accuracy.pdf
│ │ ├── prem_speed_accuracy_3plus1.pdf
│ │ ├── speed_accuracy_combined.pdf
│ │ ├── lri_earth.pdf
│ │ └── Notes_on_SU_n__probability_relations.pdf
│ └── tools/ # Tooling for the paper, not for the library
│ └── manuscript/ # Checks the prose of main.tex against the criteria beside them
│ ├── CHECKLIST.md # The audit criteria the paper is held to
│ ├── audit_section.py # The mechanical criteria, one section at a time
│ ├── check_apparatus_prose.py # Prose about the manuscript rather than about neutrinos
│ ├── check_convolution.py # Sentence shapes a reader has to unpick: abstract predicates and echoed heads
│ ├── check_dangles.py # Dangling last lines, over the typeset page
│ ├── check_informality.py # Colloquial register in a technical text
│ ├── check_llm_tells.py # Vocabulary and rhythms that read as machine-written
│ ├── check_numerals.py # Every quoted number, ranked by how badly it could be wrong
│ ├── check_oxford.py # Three-item lists missing the serial comma
│ ├── check_paragraphs.py # Paragraph length and sentence count, from the source
│ ├── check_scale_words.py # Adjectives of scale standing where a number should
│ ├── check_voice.py # Passive rate and first-person density, by section
│ └── check_weak_prose.py # Hedges, nominalizations, filler and vague reference
└── tests/ # Regression suite, run with pytest
├── conftest.py # Shared fixtures and path setup
├── bench/ # Fair-comparison benchmark pipeline
│ ├── artifacts/ # One JSON per measured cell, written by the run and replaced by the next one
│ ├── manifest.json # Pinned versions, build profiles, thread policy, capabilities
│ ├── README.md # Entry point: state of play and where each decision is recorded
│ ├── ADVERSARIAL.md # What the adversarial check must try to break
│ ├── FINDINGS.md # What that check broke, and how each claim was verified
│ ├── OBJECTIONS.md # Each objection, the measurement answering it, and its test
│ ├── requirements.lock # Exact Python versions the benchmark venv installs
│ ├── build.sh # Clones, hash-verifies and builds all seven at their pins
│ ├── bench.hpp # The C++ harness: owns main(), every clock, the statistics
│ ├── machine.py # Environment capture, the canary, and the rejection rule
│ ├── conversions.py # The one place a physical conversion factor is computed
│ ├── gen_conversions.py # Emits conversions.h so no adapter carries a physical literal
│ ├── adapters/ # One adapter per code, physics only
│ │ ├── nufast_earth.cpp # NuFast-Earth, batched over energy and zenith, invariants hoisted
│ │ ├── nufast_lbl.cpp # NuFast-LBL, one batched call over the whole energy vector
│ │ ├── prob3.cpp # Prob3++, looped in C++ because that is its interface
│ │ ├── globes.cpp # GLoBES, looped, chord decomposition hoisted
│ │ ├── nuoscprobexact.py # This library, batched
│ │ ├── second_order.py # The second-order expansion in alpha and s13, as a measured point
│ │ ├── nusquids.py # nuSQuIDS through its multiple-energy constructor
│ │ └── nucraft.py # nuCraft, whose batching is interface-only
│ ├── runner.py # The Python harness; the only place Python timing happens
│ ├── reference.py # One 50-digit reference per code, in that code's own convention
│ ├── audit_reference.py # Measures whether each reference matches the code it judges
│ ├── sweep_accuracy.py # Every code against its own reference, over every precision knob
│ ├── run_all.py # The run matrix: which code, which grid, which protocol, which knob
│ ├── check_neutrality.py # Checks that a harness change moved the spread and not the speed
│ ├── emit_figures.py # Turns the measured artifacts into the files the paper figures draw
│ ├── audit_objections.py # Checks the measured output against each objection, one at a time
│ ├── accuracy_const.json # Accuracy of every code on the constant-density grid, per knob
│ ├── accuracy_chord.json # Accuracy of every code on the Earth chord, per knob
│ ├── accuracy_discretisation.json # Accuracy against the layer count, for the code whose precision knob is something else
│ ├── earth_plane.json # The Earth speed-accuracy plane, three flavors, for the figure
│ ├── scan_sensitivity.json # What each code costs when a fit moves Dmsq31 rather than delta_CP
│ ├── speed_accuracy_plane.json # Every timed point paired with the accuracy its own setting reached
│ └── reference_audit.json # What that audit measured; the notebook renders it, never types it
├── test_bench_pipeline.py # The benchmark pipeline's fairness invariants
├── test_bench_behaviour.py # What an adapter does, not where its constants came from
├── gen_stiff_reference.py # Regenerates the fifty-digit four-flavor oracle
├── stiff_reference.json # That oracle, frozen: nine Hamiltonians in hexadecimal floats
├── test_su3_algebra.py # d tensor, star product, SU(3) invariants
├── test_oscprob4nu.py # SU(4) algebra, quartic roots, 3+1 physics
├── test_evolution_operator.py # U against an independent matrix exponential
├── test_probabilities.py # Normalization, positivity, P = |U|^2
├── test_hamiltonians.py # Sample Hamiltonians and sign conventions
├── test_reference_formulas.py # Exact result against the standard formulas
├── test_matter_eigenvalues.py # Matter spectrum, against Zaglauer-Schwarzer
├── test_edge_cases.py # Degenerate and near-degenerate Hamiltonians
├── test_docstrings.py # Runs the examples embedded in the docstrings
├── test_vectorized.py # The batched path, against the scalar one
├── test_vectorized_hamiltonians.py # Hamiltonians built for an array of energies
├── test_annotations.py # Annotations, and their agreement with the docs
├── test_fastkernels.py # Both backends, against each other
├── test_physical_scales.py # Both backends at the scales actually used
├── test_slabs.py # Slab composition, against expm
├── test_earth.py # PREM, geometry, and Earth probabilities
├── test_batching_and_tolerance.py # Energy batching, and the rtol/atol refinement
├── bit_capture.py # Exact-bit capture, for refactors meant to be invisible
├── bit_compare.py # Compares two captures, in ulps
├── test_documented_figures.py # Keeps the quoted performance figures agreeing
├── test_paper.py # Keeps the paper agreeing with the repository
├── test_version_consistency.py # Keeps the version agreeing wherever it is implied
├── test_nusquids_comparison.py # Against nuSQuIDS, an independent external code
├── nusquids_reference.py # Regenerates the frozen nuSQuIDS reference data
├── nusquids_reference.json # Those reference values, with their provenance
├── nusquids_scan.py # Regenerates the frozen nuSQuIDS energy scan
├── nusquids_scan.json # That scan, for the paper's figures
├── nufast_scan.json # NuFast-LBL at two Newton settings
├── speed_accuracy.json # The six-code constant-density speed-accuracy plane
├── timing_other_codes.json # Timings behind the performance figure
├── measure_performance_scaling.py # Freezes the three call routes against stack size, on a quiet machine
├── performance_scaling.json # That measurement, which the performance figure reads
├── prem_scan.py # Regenerates the two Earth speed-accuracy planes
├── prem_speed_accuracy.json # Those planes, at three flavors and at 3+1
├── const_density_scan.json # The exact reference for the comparison figure
├── external_drivers/ # Drivers for the codes that cannot be called from Python
│ ├── nufast_scan.cpp # Regenerates the NuFast-LBL curve of the comparison figure
│ ├── README.md # Every convention each one had to be told, and why
│ ├── gen_prem_header.py # Emits this library's PREM as a C header
│ ├── nufast_drv.cpp # NuFast-Earth on the PREM chord
│ ├── nufast_earth_prem.txt
│ ├── globes_drv.c # GLoBES on the PREM chord
│ ├── globes_prem.txt
│ ├── prob3_drv.cpp # Prob3++ on the PREM chord
│ ├── prob3_prem.txt
│ ├── perf_nusquids.py # nuSQuIDS on the constant-density scan
│ ├── globes_perf.c # GLoBES on the same scan
│ ├── prob3_perf.cpp # Prob3++ on the same scan
│ ├── nufast_lbl_perf.cpp # NuFast-LBL at two Newton settings
│ ├── const_scan.py # The 50-digit reference for the comparison figure
│ ├── globes_scan.c # GLoBES on that comparison
│ └── prob3_scan.cpp # Prob3++ on that comparison
└── test_file_tree.py # Keeps this tree in step with the repository
The three core modules are self-contained --- they need only numpy and the standard library --- so copying src/oscprob2nu.py, src/oscprob3nu.py or src/oscprob4nu.py into your own project is a supported way to use NuOscProbExact. Each imports fastkernels if it is available and does without it if it is not, so a lone copy works and simply runs the NumPy path; a test copies each of the three out and exercises it that way. Adding src/ to the path works too, and is what the bundled examples do:
import sys
sys.path.append('/path/to/NuOscProbExact/src')
import oscprob3nuRun the worked examples.
Inside the directory examples/, we provide several example files to get you started. Each is runnable as it stands and prints the probabilities it computes; Usage and examples below walks through the first of them. To run any of the examples, just execute, e.g.,
python example_2nu_trivial.pyInspecting the example files and reading their description below will help you to learn how to use NuOscProbExact in your own project.
Renamed: this directory was called
test/in version 1.0.0 of the code, and is named that way in version 2 of the paper. It becameexamples/to stop it being confused withtests/, which holds the regression suite.
Run the regression tests.
cd /path/to/NuOscProbExact
pytestThese check the SU(2), SU(3) and SU(4) machinery against independent computations --- unitarity of the evolution operator, agreement with scipy.linalg.expm, agreement with the standard oscillation formulas, and the sign conventions of the sample Hamiltonians --- and run every example embedded in the docstrings.
Open the notebooks.
cd /path/to/NuOscProbExact
pip install -e ".[notebooks]"
jupyter lab notebooks/Twenty worked notebooks, numbered in reading order, covering the probabilities against baseline and against energy, matter and new physics, oscillograms, bi-probability plots, the Earth, arbitrary matter profiles, performance, the paper's own figures, the textbook approximations, mass ordering and the octant, antineutrinos, solar neutrinos, numerical edge cases, four-neutrino 3+1 scenarios, cross-checks with other public codes, the evolution operator itself, animated parameter sweeps, and a Hamiltonian of your own carried through a varying profile. They carry their figures inline, so they can also just be read on GitHub.
The probabilities are computed from a closed form, so a single one is quick — about 8 µs for three flavors and 1 µs for two. Most real use, though, is a scan: a curve versus baseline or energy, or an oscillogram over both. Two things make those much faster, and neither changes the answers.
Every core routine accepts a stack of Hamiltonians, an array of baselines, or both, and evaluates the whole thing in one call:
# instead of this
prob = [oscprob3nu.probabilities_3nu(h_vacuum, l) for l in baselines]
# do this
prob = oscprob3nu.probabilities_3nu(h_vacuum, baselines) # (N, 9)The sample Hamiltonians take an array of energies too, so a scan in matter is two calls and no Python loop:
h_stack = hamiltonians3nu.hamiltonian_3nu_matter(h_vacuum_energy_indep,
energies, VCC_EARTH_CRUST)
prob = oscprob3nu.probabilities_3nu(h_stack, baseline)This is the single biggest win — roughly 20–90× — and it needs no extra dependency. It works because the expansion's expensive part, the characteristic equation whose roots give the oscillation phases, depends on the Hamiltonian alone: a scan over baselines solves it once rather than once per point.
Nothing to install and nothing in your code to change: Numba comes with the package as of 1.13.0. The batched paths run as compiled machine-code loops spread over your cores instead of as a chain of NumPy array operations. Where Numba is unavailable, or switched off with fastkernels.USE_NUMBA = False, the NumPy path is used and the results are the same to round-off.
Measured on 2000-point scans, against the equivalent Python loop:
| Scan | loop | arrays | arrays + Numba |
|---|---|---|---|
| Three-flavor, vs. baseline | 38 ms | 1.8 ms (~21×) | 0.31 ms (~120×) |
| Three-flavor, vs. energy | 34 ms | 1.5 ms (~23×) | 0.20 ms (~170×) |
| Three-flavor oscillogram, 100×100 | 197 ms | 5.3 ms (~37×) | 0.85 ms (~230×) |
| Two-flavor, vs. baseline | 6.9 ms | 0.07 ms (~99×) | not used — see below |
Best of seven runs, interleaved, on one machine. These are indicative, not precise: repeated runs vary by tens of per cent, so treat them as orders of magnitude. Notebook 09 measures the same comparison on whatever machine runs it, which is the number to trust.
Checking the input costs more than the arithmetic on a large scan. Every entry point verifies that the Hamiltonian is Hermitian, because one that is not returns probabilities that still sum to one — so nothing downstream reveals the mistake. Validating a stack is a pass over it, the same order of work as evaluating it: 1.3× to 1.8× at two thousand points, and 3.2× to 5.7× at two hundred thousand, where the compiled kernel has made the evaluation fast enough that the check dominates. If your Hamiltonians come from a construction you already trust — everything hamiltonians2nu, hamiltonians3nu and hamiltonians4nu build is Hermitian to round-off — decline it with oscprob3nu.CHECK_HERMITICITY = False, and likewise on the other two modules.
Layered matter gets its own kernel. slabs and earth compose evolution operators rather than computing probabilities, so until 1.12.0 they could not use the backend at all and quietly ran the NumPy path. They now compose the whole trajectory in one compiled pass — never materialising the stack, never dispatching a product per slab — which is worth ~12×, ~12× and ~9× on a single Earth crossing at two, three and four flavors. The threshold there is one slab, not MIN_BATCH: the comparison is against a Python loop, so the kernel is ahead at every length.
An Earth scan is one call. Both the energy and the zenith angle may be arrays, and they broadcast, so an oscillogram is probabilities_3nu_earth(h, energies[None, :], costhz[:, None]) rather than a loop over its points. The geometry and the matter potentials depend on the angle alone, so a scan builds them once instead of once per energy; the chord kernels then build each slab's Hamiltonian as they go, so the stack of operators that made the batched path memory-bound is never allocated. Over two thousand energies:
| Earth scan, 2000 energies | NumPy loop | Compiled loop | Compiled array |
|---|---|---|---|
| Two flavors | 542 ms | 44 ms | 1.2 ms |
| Three flavors | 895 ms | 77 ms | 7.5 ms |
| Four flavors | 2224 ms | 272 ms | 25 ms |
That is ~38×, ~10× and ~11× over the compiled loop, and ~460×, ~120× and ~87× over the NumPy one. A 100 × 100 oscillogram at three flavors takes 42 ms against 371 ms for the loop, about 4 µs a point.
A chord is a palindrome. A neutrino crossing a spherically symmetric Earth meets every radius on the way in and again on the way out, so slab j and slab n−1−j carry the same Hamiltonian, the same width, and therefore the same operator. Composing from the centre outwards computes each of them once instead of twice — worth ~1.4×, ~1.5× and ~1.8× of the figures above. Nothing about it is specific to PREM: fastkernels.palindromic decides, and slabs asks it of any sequence, so a symmetric castle wall is composed at the same discount. fastkernels.USE_PALINDROME = False asks for the plain left-to-right product instead.
Ask for an accuracy, not a slab count. n_slabs_per_segment fixes the discretisation rather than the error, and the two are not the same thing: the error is strongly energy- and angle-dependent. Pass rtol or atol to any earth entry point and the chord is refined until the measured error meets it, on every returned probability, raising rather than silently returning something coarser when it cannot. earth.slabs_for_tolerance runs the same search on its own, which is the cheap way to set one subdivision for a whole scan. The same machinery serves profiles that are not the Earth through slabs.probabilities_{2,3,4}nu_profile, which take the profile as a callable — scaled by the baseline rather than by the last sampled position, which moves as the refinement doubles.
The backend is not used where it would not help. For three flavors it wins at every stack size, by between two and sixteen times. For two flavors it does not: that expansion reduces to a square root and a sine per element, which NumPy already does about as well as compiled code can, and the kernel additionally has to materialise the Hamiltonian stack. Below fifty thousand elements the NumPy path is quicker, so it is kept; above, the kernel leads by about 1.3–1.8×. The thresholds are measured, and the library picks whichever is faster without you doing anything.
Two costs, so the trade is visible: importing Numba takes about 140 ms against 65 ms for NumPy alone, and the first call into a kernel compiles, which takes a few seconds. The kernels are cached on disk, so later runs start in milliseconds. Neither cost falls on a single scalar probability, which never enters a kernel at all — it is the closed form, and it returns in microseconds whether or not Numba is present.
That is why this was an optional extra until 1.13.0, and the argument that changed: the costs are paid once and are invisible next to a 12× Earth crossing, while an extra is paid attention to only by readers who get this far.
- Short stacks. Below thirteen elements at three flavors, and twelve at two, the array machinery costs more than it saves, so those are evaluated one at a time automatically.
- The scalar path. It is deliberately left uncompiled: 8 µs is not worth a compilation pause on a first call.
- Turning Numba off.
fastkernels.USE_NUMBA = Falseforces the NumPy path, which is how the test suite checks that the two agree.
One thing that is worth doing by hand: build the energy-independent part of the vacuum Hamiltonian once, outside any scan, since it does not depend on the energy. The bundled examples all do this.
There are three core modules, one per flavor count: oscprob2nu.py, oscprob3nu.py and oscprob4nu.py. Each is stand-alone apart from the dependencies described above. Install the package, or add src/ to the path, which is what the bundled examples do:
import sys
sys.path.append('../src')What follows is the short version: what the functions take and return, and four examples that between them cover a single probability, a whole scan, four flavors, and your own Hamiltonian. Everything else --- matter, non-standard interactions, Lorentz-invariance violation, oscillograms, the Earth, the expansion coefficients --- lives in the runnable scripts in examples/ and in the notebooks, which store their figures inline and are executed by CI. It is not repeated here, so there is one copy of each to keep correct.
Most of the time you want probabilities, not the intermediate steps. The routine is probabilities_Nnu in oscprobNnu, and it takes a Hermitian matrix and a baseline:
| Flavors | Call | Returns |
|---|---|---|
| 2 | oscprob2nu.probabilities_2nu(h, L) |
4 values: Pee, Pem, Pme, Pmm |
| 3 | oscprob3nu.probabilities_3nu(h, L) |
9 values: Pee, Pem, Pet, Pme, ..., Ptt |
| 4 | oscprob4nu.probabilities_4nu(h, L) |
16 values: Pee, Pem, Pet, Pes, Pme, ..., Pss |
In every case the initial flavor varies slowest, so P[n*alpha + beta] is P(nu_alpha -> nu_beta) for n flavors. The two-flavor labels could equally be Pmm, Pmt, Ptm, Ptt --- which pair of flavors they describe is set by the Hamiltonian you pass, not by the code.
The evolution operator itself is available too, as evolution_operator_Nnu(h, L), if you want to compose it across segments or propagate a density matrix rather than read off probabilities.
Important: The Hamiltonian must be Hermitian, and every entry point checks that it is: one that is not raises
ValueErrorrather than returning numbers. The check is there because the numbers it would otherwise return still sum to one, so nothing downstream would reveal the mistake. Performance gives what the check costs and how to decline it where your Hamiltonians are Hermitian by construction.
About the units: These modules assume no units for any of the model parameters, so you need to pass values with consistent ones --- all that is required is that
H*Lbe dimensionless. The moduleglobaldefsprovides physical constants and conversion factors, includingCONV_KM_TO_INV_EV, which converts a baseline in km to eV^{-1}.
Three-flavor oscillations in vacuum, at a fixed energy and baseline. hamiltonian_3nu_vacuum_energy_independent returns the vacuum Hamiltonian without the 1/E prefactor, so that a scan over energies computes it once and divides by a varying E:
# Find this example in NuOscProbExact/examples/example_3nu_vacuum.py
import numpy as np
import oscprob3nu
import hamiltonians3nu
from globaldefs import *
energy = 1.e9 # Neutrino energy [eV]
baseline = 1.3e3 # Baseline [km]
# NuFit best-fit mixing parameters, pulled from globaldefs. NO means
# "normal ordering"; change NO to IO for inverted ordering
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.asarray(h_vacuum_energy_indep)/energy
Pee, Pem, Pet, Pme, Pmm, Pmt, Pte, Ptm, Ptt = oscprob3nu.probabilities_3nu(
h_vacuum, baseline*CONV_KM_TO_INV_EV)
print("Pee = %6.5f, Pem = %6.5f, Pet = %6.5f" % (Pee, Pem, Pet))
print("Pme = %6.5f, Pmm = %6.5f, Pmt = %6.5f" % (Pme, Pmm, Pmt))
print("Pte = %6.5f, Ptm = %6.5f, Ptt = %6.5f" % (Pte, Ptm, Ptt))This returns
Pee = 0.92768, Pem = 0.01432, Pet = 0.05800
Pme = 0.04023, Pmm = 0.37887, Pmt = 0.58090
Pte = 0.03210, Ptm = 0.60680, Ptt = 0.36110Each row sums to one, as it must.
Antineutrinos: pass
-dCPinstead ofdCP, and flip the sign of the matter potential. Notebook 13 works through both, and the two ways to get it wrong.
Do not call the routine in a Python loop. Every core routine accepts a stack of Hamiltonians, an array of baselines, or both broadcast against each other, and evaluates the lot in one call --- which is tens of times faster and gives identical results:
import numpy as np
import oscprob3nu
import hamiltonians3nu
from globaldefs import *
baseline = 1.3e3*CONV_KM_TO_INV_EV # [eV^{-1}]
energies = np.logspace(-1.0, 1.0, 200)*1.e9 # 0.1 to 10 GeV [eV]
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)
# One Hamiltonian per energy, stacked along a leading axis
h_stack = np.asarray(h_vacuum_energy_indep)/energies[:, None, None]
prob = oscprob3nu.probabilities_3nu(h_stack, baseline) # shape (200, 9)
prob_ee, prob_em, prob_et = prob[:, 0], prob[:, 1], prob[:, 2]
print("prob.shape =", prob.shape)
print("P_ee at %5.2f GeV = %.5f" % (energies[0]/1.e9, prob_ee[0]))
print("P_ee at %5.2f GeV = %.5f" % (energies[-1]/1.e9, prob_ee[-1]))This returns
prob.shape = (200, 9)
P_ee at 0.10 GeV = 0.24693
P_ee at 10.00 GeV = 0.98582The same works for a scan over baselines, or for both at once to build an oscillogram. The sample Hamiltonians in matter accept an array of energies directly, so a matter scan is two calls and no loop. Notebook 02 plots these curves, notebook 04 builds an oscillogram, and notebook 09 measures what the broadcasting buys on whatever machine runs it.
oscprob4nu works exactly the same way, with a 4x4 Hamiltonian and sixteen probabilities. With the fourth state read as sterile, the flavor order is (nu_e, nu_mu, nu_tau, nu_s):
import numpy as np
import oscprob4nu
import hamiltonians4nu
from globaldefs import *
# Three extra mixing angles and one extra mass-squared splitting,
# here Dm41^2 = 1 eV^2
h_vacuum_energy_indep = hamiltonians4nu.hamiltonian_4nu_vacuum_energy_independent(
S12_NO_BF, S23_NO_BF, S13_NO_BF,
np.sqrt(0.10), np.sqrt(0.10), 0.0,
DCP_NO_BF, D21_NO_BF, D31_NO_BF, 1.0)
prob = oscprob4nu.probabilities_4nu(
np.asarray(h_vacuum_energy_indep)/1.e9, 1.3e3*CONV_KM_TO_INV_EV)
print("%d probabilities" % len(prob))
print("Pee = %6.5f, Pem = %6.5f, Pet = %6.5f, Pes = %6.5f" % tuple(prob[0:4]))
print("they sum to %.5f" % sum(prob[0:4]))16 probabilities
Pee = 0.76700, Pem = 0.00149, Pet = 0.05220, Pes = 0.17931
they sum to 1.00000Notebook 16 works a 3+1 scenario through properly --- the sterile entry in the matter potential, a short-baseline scan, the sterile matter resonance through the Earth --- and explains why four flavors is where the closed form ends.
Nothing above is a special case in the code: vacuum, matter, non-standard interactions and Lorentz-invariance violation are each just a different Hermitian matrix handed to the same routine. So your own model is too. Usually you will want to add a term to the vacuum Hamiltonian, where hamiltonian_mymodel is yours to write and returns a 3x3 matrix:
import numpy as np
import oscprob3nu
import hamiltonians3nu
from globaldefs import *
energy = 1.e9 # Neutrino energy [eV]
baseline = 1.3e3 # Baseline [km]
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.asarray(h_vacuum_energy_indep)/energy
h_mymodel = h_vacuum + hamiltonian_mymodel(mymodel_parameters)
prob = oscprob3nu.probabilities_3nu(h_mymodel, baseline*CONV_KM_TO_INV_EV)hamiltonian_mymodel may depend on the energy too. For two flavors it returns a 2x2 matrix instead, and for four, a 4x4 one. Passing an arbitrary matrix directly, with no vacuum term at all, works exactly as you would expect --- see examples/example_3nu_trivial.py.
Each of these is a runnable script; none of them is transcribed into this file, so there is a single copy to keep correct.
| Scenario | Script | Notebook |
|---|---|---|
| Arbitrary Hamiltonian, 2 and 3 flavors | example_2nu_trivial.py, example_3nu_trivial.py |
01 |
| Vacuum, 2 and 3 flavors | example_2nu_vacuum.py, example_3nu_vacuum.py |
02 |
| Constant-density matter | example_3nu_matter.py |
03 |
| Matter with non-standard interactions | example_3nu_nsi.py |
03 |
| Lorentz-invariance violation | example_3nu_liv.py |
03 |
| SU(2) and SU(3) expansion coefficients, and the evolution operator | example_2nu_vacuum_coeffs.py, example_3nu_vacuum_coeffs.py |
18 |
| Layered matter, and the Earth through PREM | — | 06, 07 |
The numerical recipes page collects the same material as runnable snippets, the tutorial notebooks page indexes the twenty in reading order, and the API reference documents every routine, with examples that are executed when the documentation is built rather than pasted beside it.
Twenty worked notebooks live in notebooks/, numbered in reading order. They carry their figures inline, so they render on GitHub without being run. The documentation lists them with the same descriptions, grouped by what they are for, on its tutorial notebooks page:
| Notebook | What it covers |
|---|---|
| 01 Basics | Units, one probability, and why to pass arrays rather than loop |
| 02 Vacuum oscillations | Against baseline and against energy, checked against the textbook formula |
| 03 Matter, NSI, LIV | Constant-density matter and two kinds of new physics |
| 04 Oscillograms | A 240x240 energy-baseline map in a single call |
| 05 Bi-probability | CP ellipses, in vacuum and in matter |
| 06 The Earth and PREM | The density profile, chord geometry, slabs, and their convergence |
| 07 Through the Earth | Zenith-angle scans, an Earth oscillogram, and real baselines |
| 08 Unusual density profiles | Castle-wall and serrated profiles, and parametric enhancement |
| 09 Performance | Looping versus broadcasting, and the compiled backend, measured live |
| 10 The paper's figures | The two figures from arXiv:1904.12391 |
| 11 Exact vs approximations | Where the familiar formulas agree, and where they do not |
| 12 Ordering and octant | Normal against inverted, and the θ₂₃ octant degeneracy |
| 13 Antineutrinos | Conjugate and flip the potential — and two ways to get it wrong |
| 14 Solar and the MSW resonance | The adiabatic resonance, validated — and why slabs are the wrong tool for it |
| 15 Numerical edge cases | Degenerate spectra, and what returns a number instead of NaN |
| 16 Four neutrinos | A 3+1 sterile state through SU(4), and why the method stops at four |
| 17 Cross-checks with other codes | Corroboration from nuSQuIDS and from Zaglauer–Schwarzer, and the conventions that had to be matched |
| 18 The evolution operator | The operator itself, the group property, and the SU(n) coefficients underneath |
| 19 Animated scenes | Four parameter sweeps drawn as stills, and the reel they came from |
| 20 An arbitrary Hamiltonian | A Hamiltonian of your own carried through three varying profiles, and how to put one through the Earth |
Run them with pip install -e ".[notebooks]" and jupyter lab notebooks/. Every one of them is executed by CI, so an example that stops working fails the build.
All of the modules provided in NuOscProbExact have been documented using Python docstrings, written in numpydoc format so that they can be rendered directly by Sphinx with the numpydoc extension. They are human-readable by opening the source .py files. Alternatively, they can be printed from within an interactive Python session.
Every Examples block in the docstrings is executed when the documentation is built, so the results shown on the API page are produced by the code rather than pasted beside it, and cannot drift. The regression suite runs the same blocks on every supported Python (tests/test_docstrings.py), which the documentation build --- one job, one interpreter --- would not catch.
A full Sphinx project lives in docs/. Build it with
pip install -r docs/requirements.txt
cd docs && make htmland open docs/build/html/index.html. It contains an installation guide, a quickstart, a description of the method and its sign conventions, the API reference generated from the docstrings, a bibliography, and the changelog.
Notable changes between versions are recorded in CHANGELOG.md, which the documentation renders as its own page, so there is a single source of truth.
To view the documentation of a module from within an interactive Python session, run, e.g.,
import oscprob3nu
print(oscprob3nu.__doc__)This will print to screen a description of what the module does (in this example, oscprob3nu) and a list of the functions that it contains, including a description of each.
To view the documentation of a particular function from within an interactive Python session, run, e.g.,
import oscprob3nu
help(oscprob3nu.hamiltonian_3nu_coefficients)This will print to screen a description of what the function does (in the example above, oscprob3nu.hamiltonian_3nu_coefficients), a list and description of its input parameters, and a description of the values that it returns.
If you use NuOscProbExact in your work, we ask you that you please cite the following paper: Mauricio Bustamante, NuOscProbExact: a general-purpose code to compute exact two-flavor and three-flavor neutrino oscillation probabilities (arXiv:1904.12391). The papers that cite it are listed on INSPIRE.
If you are citing NuOscProbExact in a document that will be uploaded to the arXiv, please consider using the LaTeX or BibTeX entries provided by INSPIRE (link here):
@article{Bustamante:2019ggq,
author = "Bustamante, Mauricio",
title = "{NuOscProbExact: a general-purpose code to compute
exact two-flavor and three-flavor neutrino
oscillation probabilities}",
year = "2019",
eprint = "1904.12391",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
SLACcitation = "%%CITATION = ARXIV:1904.12391;%%"
}
The paper covers two and three flavors, which is what it was written about. To cite the software — a particular version of it, including the four-flavor extension that came after the paper — use the Zenodo DOI badge at the top of this file, which resolves to the most recently archived release. Zenodo mints a DOI per GitHub Release, so citing a specific version means citing the DOI archived for it.
NuOscProbExact is released under the MIT License. The full text ships with the source, as LICENSE in the repository root.