Metadata-Version: 2.4
Name: nuoscprobexact
Version: 1.11.0
Summary: Exact two-, three- and four-flavor neutrino oscillation probabilities for arbitrary time-independent Hamiltonians
Author-email: Mauricio Bustamante <mbustamante@gmail.com>
License: MIT License
        
        Copyright (c) 2019 Mauricio Bustamante
        
        Permission is hereby granted, free of charge, to any person obtaining a copy
        of this software and associated documentation files (the "Software"), to deal
        in the Software without restriction, including without limitation the rights
        to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
        copies of the Software, and to permit persons to whom the Software is
        furnished to do so, subject to the following conditions:
        
        The above copyright notice and this permission notice shall be included in all
        copies or substantial portions of the Software.
        
        THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
        IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
        FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
        AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
        LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
        OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
        SOFTWARE.
        
Project-URL: Homepage, https://github.com/mbustama/NuOscProbExact
Project-URL: Documentation, https://mbustama.github.io/NuOscProbExact/
Project-URL: Changelog, https://github.com/mbustama/NuOscProbExact/blob/main/CHANGELOG.md
Project-URL: Issues, https://github.com/mbustama/NuOscProbExact/issues
Project-URL: Paper, https://arxiv.org/abs/1904.12391
Keywords: neutrino,oscillations,physics,sterile neutrinos,SU(2),SU(3),SU(4)
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy
Provides-Extra: fast
Requires-Dist: numba; extra == "fast"
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Requires-Dist: scipy; extra == "test"
Requires-Dist: coverage; extra == "test"
Requires-Dist: pytest-cov; extra == "test"
Provides-Extra: docs
Requires-Dist: sphinx; extra == "docs"
Requires-Dist: sphinx-rtd-theme; extra == "docs"
Requires-Dist: sphinx-copybutton; extra == "docs"
Requires-Dist: sphinxcontrib-bibtex; extra == "docs"
Requires-Dist: numpydoc; extra == "docs"
Requires-Dist: myst-parser; extra == "docs"
Requires-Dist: jupyter-sphinx; extra == "docs"
Requires-Dist: ipykernel; extra == "docs"
Requires-Dist: matplotlib; extra == "docs"
Requires-Dist: scipy; extra == "docs"
Provides-Extra: notebooks
Requires-Dist: jupyter; extra == "notebooks"
Requires-Dist: nbclient; extra == "notebooks"
Requires-Dist: nbformat; extra == "notebooks"
Requires-Dist: ipykernel; extra == "notebooks"
Requires-Dist: matplotlib; extra == "notebooks"
Requires-Dist: scipy; extra == "notebooks"
Dynamic: license-file

[![tests](https://github.com/mbustama/NuOscProbExact/actions/workflows/tests.yml/badge.svg)](https://github.com/mbustama/NuOscProbExact/actions/workflows/tests.yml)
[![Code Quality](https://github.com/mbustama/NuOscProbExact/actions/workflows/lint.yml/badge.svg)](https://github.com/mbustama/NuOscProbExact/actions/workflows/lint.yml)
[![codecov](https://codecov.io/gh/mbustama/NuOscProbExact/branch/main/graph/badge.svg)](https://codecov.io/gh/mbustama/NuOscProbExact)
[![Documentation](https://img.shields.io/badge/docs-GitHub%20Pages-blue.svg)](https://mbustama.github.io/NuOscProbExact/)
[![PyPI](https://img.shields.io/pypi/v/nuoscprobexact.svg)](https://pypi.org/project/nuoscprobexact/)
[![Downloads](https://pepy.tech/badge/nuoscprobexact)](https://pepy.tech/project/nuoscprobexact)
[![arXiv](https://img.shields.io/badge/arXiv-1904.12391-orange.svg)](https://arxiv.org/abs/1904.12391)
[![DOI](https://zenodo.org/badge/182178323.svg)](https://zenodo.org/badge/latestdoi/182178323)
[![License: MIT](https://img.shields.io/badge/License-MIT-blue.svg)](https://opensource.org/licenses/MIT)
[![Python 3.9+](https://img.shields.io/badge/python-3.9+-blue.svg)](https://www.python.org/downloads/)
[![Code style: ruff](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/astral-sh/ruff/main/assets/badge/v2.json)](https://github.com/astral-sh/ruff)

# NuOscProbExact
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 with `pytest`, and the results are cross-checked against [nuSQuIDS](https://github.com/arguelles/nuSQuIDS), an independent external code — see [notebook 17](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/17_cross_checks.ipynb).

## What you can compute

Every figure below is produced by a notebook in [`notebooks/`](https://github.com/mbustama/NuOscProbExact/tree/main/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](https://mbustama.github.io/NuOscProbExact/recipes.html) page.

| | |
|:--:|:--:|
| <img src="https://raw.githubusercontent.com/mbustama/NuOscProbExact/main/img/gallery/gallery_vacuum.png" width="380"/><br/>**Oscillation probabilities** against baseline or energy, for two, three or four flavors.<br/>[notebook 02](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/02_vacuum_oscillations.ipynb) | <img src="https://raw.githubusercontent.com/mbustama/NuOscProbExact/main/img/gallery/gallery_matter.png" width="380"/><br/>**Matter, NSI and Lorentz-invariance violation** — each just a different Hermitian matrix.<br/>[notebook 03](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/03_matter_nsi_liv.ipynb) |
| <img src="https://raw.githubusercontent.com/mbustama/NuOscProbExact/main/img/gallery/gallery_oscillogram.png" width="380"/><br/>**Oscillograms** over energy and baseline: 57 600 probabilities in a single call.<br/>[notebook 04](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/04_oscillogram.ipynb) | <img src="https://raw.githubusercontent.com/mbustama/NuOscProbExact/main/img/gallery/gallery_biprobability.png" width="380"/><br/>**CP violation**, as bi-probability ellipses in vacuum and in matter.<br/>[notebook 05](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/05_biprobability.ipynb) |
| <img src="https://raw.githubusercontent.com/mbustama/NuOscProbExact/main/img/gallery/gallery_prem.png" width="380"/><br/>**The Earth's density**, from the Preliminary Reference Earth Model.<br/>[notebook 06](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/06_earth_and_prem.ipynb) | <img src="https://raw.githubusercontent.com/mbustama/NuOscProbExact/main/img/gallery/gallery_earth.png" width="380"/><br/>**Neutrinos through the Earth**, in energy and zenith angle, or between two named sites.<br/>[notebook 07](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/07_earth_probabilities.ipynb) |
| <img src="https://raw.githubusercontent.com/mbustama/NuOscProbExact/main/img/gallery/gallery_profiles.png" width="380"/><br/>**Arbitrary matter profiles** — castle walls and worse, exactly.<br/>[notebook 08](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/08_unusual_density_profiles.ipynb) | <img src="https://raw.githubusercontent.com/mbustama/NuOscProbExact/main/img/gallery/gallery_ordering.png" width="380"/><br/>**Mass ordering and the θ₂₃ octant**, separated by matter through the Earth.<br/>[notebook 12](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/12_ordering_and_octant.ipynb) |
| <img src="https://raw.githubusercontent.com/mbustama/NuOscProbExact/main/img/gallery/gallery_sterile.png" width="380"/><br/>**Four flavors: a 3+1 sterile state**, resolved at a short baseline.<br/>[notebook 16](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/16_four_neutrinos.ipynb) | <img src="https://raw.githubusercontent.com/mbustama/NuOscProbExact/main/img/gallery/gallery_sterile_earth.png" width="380"/><br/>**The sterile matter resonance** through the Earth, in energy and zenith angle.<br/>[notebook 16](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/16_four_neutrinos.ipynb) |

## Contents

1. [What you can compute](#what-you-can-compute)

2. [What is NuOscProbExact?](#what-is-nuoscprobexact)
   1. [What it does](#what-it-does)
   2. [What it does not do](#what-it-does-not-do)
   3. [When to use Magnus instead](#when-to-use-magnus-instead)

3. [Requirements](#requirements)

4. [Installation](#installation)

5. [Performance](#performance)

6. [Usage and examples](#usage-and-examples)
   1. [Basics](#basics)
   2. [A first probability](#a-first-probability)
   3. [Whole scans in one call](#whole-scans-in-one-call)
   4. [Four flavors: a 3+1 sterile state](#four-flavors-a-31-sterile-state)
   5. [Arbitrary Hamiltonians](#arbitrary-hamiltonians)
   6. [Where the rest is](#where-the-rest-is)

7. [Notebooks](#notebooks)

8. [Documentation and help](#documentation-and-help)

9. [Citing](#citing)

10. [License](#license)


## What is NuOscProbExact?

**NuOscProbExact** is a Python implementation of the method developed by [Ohlsson & Snellman](https://arxiv.org/abs/hep-ph/9910546) 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](https://arxiv.org/abs/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](https://mbustama.github.io/NuOscProbExact/methodology.html#why-the-method-stops-at-four).

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](https://mbustama.github.io/NuOscProbExact/methodology.html).

### What it does

* **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](https://github.com/arguelles/nuSQuIDS) to round-off once conventions are matched ([notebook 17](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/17_cross_checks.ipynb)).  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.**  `oscprob4nu` carries 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.**  `slabs` propagates 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.**  `earth` builds 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 and `earth` builds it per slab.
* **An optional compiled backend.**  With `numba` installed, large batched calls run on compiled kernels; without it the NumPy path is used and the answers are the same to round-off.

### What it does not do

* **Hamiltonians that vary continuously along the trajectory.**  See [When to use Magnus instead](#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.

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

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 ~10<sup>4</sup> steps per resonance crossing |
| Open systems: decay, decoherence | neither | Needs a Lindblad solver, not a unitary one |

[Notebook 14](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/14_solar_and_adiabatic_msw.ipynb) 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](#citing).


## Requirements

**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.  Where a row names an extra, that is the one to install; the rest need nothing beyond `numpy`.  The commands are under [Installation](#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 *(optional)* | `numba` | `fast` |
| 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` |

Only `numpy` is ever required.  `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` is entirely optional — it is worth roughly 1.5x to 20x on large scans, depending on their size and the number of flavors, and without it the NumPy path is used and the results are identical to round-off.


## Installation

**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 optional `numba` backend still has a wheel.

### From PyPI (recommended)

```shell
pip install nuoscprobexact
```

That is the whole installation.  The only required dependency is `numpy`.

The optional extras add what each task needs, and can be combined:

```shell
pip install "nuoscprobexact[fast]"       # numba, for the compiled batched kernels
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 documentation
```

Then, in your own code:

```python
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.

### From GitHub

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:

```shell
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:

```text
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
├── LICENSE                          # MIT license
├── README.md                        # Project overview and worked examples
├── pyproject.toml                   # Packaging metadata and pytest configuration
├── 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
│       ├── 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
│   ├── 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_biprobability.png
│       ├── gallery_earth.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
│   └── make_notebooks.py            # Generates and executes all of the above
├── 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               # Optional Numba kernels, with a NumPy fallback
│   ├── slabs.py                     # Propagation across adjacent slabs
│   └── earth.py                     # PREM, chord geometry, and Earth crossings
└── tests/                           # Regression suite, run with pytest
    ├── conftest.py                  # Shared fixtures and path setup
    ├── 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_documented_figures.py   # Keeps the quoted performance figures agreeing
    ├── 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
    └── test_file_tree.py            # Keeps this tree in step with the repository
```

### Without installing anything

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:

```python
import sys
sys.path.append('/path/to/NuOscProbExact/src')

import oscprob3nu
```

### Checking the installation

**Run 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](#usage-and-examples) below walks through the first of them.  To run any of the examples, just execute, *e.g.*,
   ```shell
   python example_2nu_trivial.py
   ```
   Inspecting 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](https://arxiv.org/abs/1904.12391).  It became `examples/` to stop it being confused with `tests/`, which holds the regression suite.

**Run the regression tests.**
   ```shell
   cd /path/to/NuOscProbExact
   pytest
   ```
   These 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.**
   ```shell
   cd /path/to/NuOscProbExact
   pip install -e ".[notebooks]"
   jupyter lab notebooks/
   ```
   Eighteen 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, and the evolution operator itself.  They carry their figures inline, so they can also just be read on GitHub.

## Performance

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.

### 1. Pass arrays instead of looping

Every core routine accepts a **stack** of Hamiltonians, an array of baselines, or both, and evaluates the whole thing in one call:

```python
# 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:

```python
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.

### 2. Install Numba, if the scans are large

```shell
pip install "nuoscprobexact[fast]"
```

Nothing in your code changes.  If [Numba](https://numba.pydata.org) is importable, the batched paths run as compiled machine-code loops spread over your cores instead of as a chain of NumPy array operations; if it is not, 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](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/09_performance.ipynb) 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.

**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 compiles, which takes a few seconds.  The kernels are cached on disk, so later runs start in milliseconds.  This is why it is an optional extra and not a dependency.

### What you do not have to think about

* **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 = False` forces 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.

## Usage and examples

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](#requirements).  Install the package, or add `src/` to the path, which is what the bundled examples do:

```python
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/`](https://github.com/mbustama/NuOscProbExact/tree/main/examples/) and in the [notebooks](#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.


### Basics

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 `ValueError` rather 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](#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*L` be dimensionless.  The module `globaldefs` provides physical constants and conversion factors, including `CONV_KM_TO_INV_EV`, which converts a baseline in km to eV^{-1}.


### A first probability

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*:

```python
# 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

```shell
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.36110
```

Each row sums to one, as it must.

> **Antineutrinos:** pass `-dCP` instead of `dCP`, and flip the sign of the matter potential.  [Notebook 13](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/13_antineutrinos.ipynb) works through both, and the two ways to get it wrong.


### Whole scans in one call

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](#performance) and gives identical results:

```python
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

```shell
prob.shape = (200, 9)
P_ee at  0.10 GeV = 0.24693
P_ee at 10.00 GeV = 0.98582
```

The 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](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/02_vacuum_oscillations.ipynb) plots these curves, [notebook 04](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/04_oscillogram.ipynb) builds an oscillogram, and [notebook 09](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/09_performance.ipynb) measures what the broadcasting buys on whatever machine runs it.


### Four flavors: a 3+1 sterile state

`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):

```python
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]))
```

```shell
16 probabilities
Pee = 0.76700, Pem = 0.00149, Pet = 0.05220, Pes = 0.17931
they sum to 1.00000
```

[Notebook 16](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/16_four_neutrinos.ipynb) 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.


### Arbitrary Hamiltonians

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:

```python
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`](https://github.com/mbustama/NuOscProbExact/blob/main/examples/example_3nu_trivial.py).


### Where the rest is

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`](https://github.com/mbustama/NuOscProbExact/blob/main/examples/example_2nu_trivial.py), [`example_3nu_trivial.py`](https://github.com/mbustama/NuOscProbExact/blob/main/examples/example_3nu_trivial.py) | [01](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/01_basics.ipynb) |
| Vacuum, 2 and 3 flavors | [`example_2nu_vacuum.py`](https://github.com/mbustama/NuOscProbExact/blob/main/examples/example_2nu_vacuum.py), [`example_3nu_vacuum.py`](https://github.com/mbustama/NuOscProbExact/blob/main/examples/example_3nu_vacuum.py) | [02](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/02_vacuum_oscillations.ipynb) |
| Constant-density matter | [`example_3nu_matter.py`](https://github.com/mbustama/NuOscProbExact/blob/main/examples/example_3nu_matter.py) | [03](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/03_matter_nsi_liv.ipynb) |
| Matter with non-standard interactions | [`example_3nu_nsi.py`](https://github.com/mbustama/NuOscProbExact/blob/main/examples/example_3nu_nsi.py) | [03](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/03_matter_nsi_liv.ipynb) |
| Lorentz-invariance violation | [`example_3nu_liv.py`](https://github.com/mbustama/NuOscProbExact/blob/main/examples/example_3nu_liv.py) | [03](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/03_matter_nsi_liv.ipynb) |
| SU(2) and SU(3) expansion coefficients, and the evolution operator | [`example_2nu_vacuum_coeffs.py`](https://github.com/mbustama/NuOscProbExact/blob/main/examples/example_2nu_vacuum_coeffs.py), [`example_3nu_vacuum_coeffs.py`](https://github.com/mbustama/NuOscProbExact/blob/main/examples/example_3nu_vacuum_coeffs.py) | [18](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/18_evolution_operator.ipynb) |
| Layered matter, and the Earth through PREM | — | [06](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/06_earth_and_prem.ipynb), [07](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/07_earth_probabilities.ipynb) |

The [numerical recipes](https://mbustama.github.io/NuOscProbExact/recipes.html) page collects the same material as runnable snippets, and the [API reference](https://mbustama.github.io/NuOscProbExact/functions.html) documents every routine, with examples that are executed when the documentation is built rather than pasted beside it.


## Notebooks

Eighteen worked notebooks live in [`notebooks/`](https://github.com/mbustama/NuOscProbExact/tree/main/notebooks/), numbered in reading order.  They carry their figures inline, so they render on GitHub without being run:

| Notebook | What it covers |
|---|---|
| [01 Basics](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/01_basics.ipynb) | Units, one probability, and why to pass arrays rather than loop |
| [02 Vacuum oscillations](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/02_vacuum_oscillations.ipynb) | Against baseline and against energy, checked against the textbook formula |
| [03 Matter, NSI, LIV](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/03_matter_nsi_liv.ipynb) | Constant-density matter and two kinds of new physics |
| [04 Oscillograms](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/04_oscillogram.ipynb) | A 240x240 energy-baseline map in a single call |
| [05 Bi-probability](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/05_biprobability.ipynb) | CP ellipses, in vacuum and in matter |
| [06 The Earth and PREM](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/06_earth_and_prem.ipynb) | The density profile, chord geometry, slabs, and their convergence |
| [07 Through the Earth](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/07_earth_probabilities.ipynb) | Zenith-angle scans, an Earth oscillogram, and real baselines |
| [08 Unusual density profiles](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/08_unusual_density_profiles.ipynb) | Castle-wall and serrated profiles, and parametric enhancement |
| [09 Performance](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/09_performance.ipynb) | Looping versus broadcasting, and the compiled backend, measured live |
| [10 The paper's figures](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/10_paper_figures.ipynb) | The two figures from [arXiv:1904.12391](https://arxiv.org/abs/1904.12391) |
| [11 Exact vs approximations](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/11_exact_vs_approximations.ipynb) | Where the familiar formulas agree, and where they do not |
| [12 Ordering and octant](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/12_ordering_and_octant.ipynb) | Normal against inverted, and the θ₂₃ octant degeneracy |
| [13 Antineutrinos](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/13_antineutrinos.ipynb) | Conjugate *and* flip the potential — and two ways to get it wrong |
| [14 Solar and the MSW resonance](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/14_solar_and_adiabatic_msw.ipynb) | The adiabatic resonance, validated — and why slabs are the wrong tool for it |
| [15 Numerical edge cases](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/15_numerical_edge_cases.ipynb) | Degenerate spectra, and what returns a number instead of NaN |
| [16 Four neutrinos](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/16_four_neutrinos.ipynb) | A 3+1 sterile state through SU(4), and why the method stops at four |
| [17 Cross-checks with other codes](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/17_cross_checks.ipynb) | Corroboration from nuSQuIDS and from Zaglauer–Schwarzer, and the conventions that had to be matched |
| [18 The evolution operator](https://github.com/mbustama/NuOscProbExact/blob/main/notebooks/18_evolution_operator.ipynb) | The operator itself, the group property, and the SU(*n*) coefficients underneath |

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.


## Documentation and help

All of the modules provided in **NuOscProbExact** have been documented using Python docstrings, written in [numpydoc](https://numpydoc.readthedocs.io/) format so that they can be rendered directly by [Sphinx](https://www.sphinx-doc.org/) 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](https://mbustama.github.io/NuOscProbExact/functions.html) 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
```shell
pip install -r docs/requirements.txt
cd docs && make html
```
and open `docs/build/html/index.html`.  It contains an installation guide, a quickstart, a description of the [method](https://github.com/mbustama/NuOscProbExact/blob/main/docs/source/methodology.rst) 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](https://github.com/mbustama/NuOscProbExact/blob/main/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.*,
```python
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.*,
```python
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.


## Citing

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](https://arxiv.org/abs/1904.12391)).

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](https://inspirehep.net/literature/1731803)):
```
@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.


## License

**NuOscProbExact** is released under the [MIT License](https://opensource.org/licenses/MIT).  The full text ships with the source, as [`LICENSE`](https://github.com/mbustama/NuOscProbExact/blob/main/LICENSE) in the repository root.

