I have posted a substantially expanded version of the NuOscProbExact paper, first made public in 2019. The original described a method for computing exact two- and three-flavor oscillation probabilities for an arbitrary time-independent Hamiltonian, without ever diagonalizing it. The paper has since roughly doubled in length, and most of what is new was not in the first version at all.

Four flavors. The method now runs at four flavors as well as two and three, so a sterile neutrino is a setting rather than a different program. The two- and three-flavor cases expand the evolution operator in Pauli and Gell-Mann matrices; four comes instead from Sylvester’s formula, the spectral form of the evolution operator, which holds at any number of flavors. Four is the largest case for which the eigenvalues can be written in radicals.

Matter of varying density, including the Earth. A density profile is cut into slabs of constant density, each solved exactly, and the operators composed in order, so the approximation is in the profile and never in the propagation through it. The code carries PREM, cuts the chord at PREM’s own density jumps, and reads the density at each slab’s midpoint, so the error falls as the square of the slab width. You state a tolerance instead of guessing a slab count.

A measured comparison against six public codes. New sections set the method beside GLoBES, Prob3++, nuCraft, NuFast-LBL, NuFast-Earth, and nuSQuIDS on two speed-against-accuracy planes, one at constant density and one along a chord through the Earth. Most of the work was in the matching: every matter potential, length convention, and oscillation parameter had to be made common first. At constant density, the probability is limited by round-off alone, four decades finer than the nearest of the others. Through the Earth, it is the quickest of the six at every accuracy but the coarsest.

What the accuracy is worth. A new section converts the accuracy axis into something physical: at 1300 km through 3 g / cm³, an error of 10⁻³ in the probability is the size of a shift of the CP-violation phase by about 2°, while 10⁻⁸ is the size of a shift of 2 × 10⁻⁵ degrees. The case for an exact probability is not that a fit needs fifteen digit; it is that an exact one carries no calculation error to bound, at any baseline, density, or set of parameters.

There are also new appendices on the structure of the code, the test suite, tutorial notebooks, and the four-flavor Hamiltonians.

On the code side—the rebuilt API, the documentation site, and the notebooks—see the earlier post below. Install NuOscProbExact via PyPI:

pip install nuoscprobexact