Methodology
This page sketches how the closed-form probabilities are obtained, and records the conventions the code uses. For the derivation, see [Bus19] and the original treatment in [OS00].
The idea
For a time-independent Hamiltonian, the evolution operator is the matrix exponential
Rather than evaluating that exponential numerically, expand \(H\) in the basis of generators of SU(2) or SU(3) and exponentiate analytically, using the closure of the algebra. The result is a closed form for \(U(L)\) in the same basis, whose coefficients are elementary functions of two invariants of the Hamiltonian.
Two flavors
Expand the Hamiltonian in the Pauli matrices,
with \(h_k\) real for Hermitian \(H\). The term \(h_0\) multiplies the identity, so it contributes only an overall phase and is dropped; the code works with the traceless part throughout, which leaves the probabilities unchanged.
Exponentiating gives
with \(|h|^2 = h_1^2 + h_2^2 + h_3^2\). The transition probability follows directly:
Both \(h_1\) and \(h_2\) appear. \(h_2 = -\mathrm{Im}(H_{12})\) vanishes only when the off-diagonal entry is real, which is why the \(h_2\) term is easy to lose and hard to notice: it makes no difference in vacuum or in matter of constant density, and all the difference for a CP-violating Hamiltonian.
When \(|h| = 0\) the Hamiltonian is proportional to the identity, there is no flavor evolution, and the code takes the limit \(\sin(|h|L)/|h| \to L\) explicitly.
Three flavors
Expand in the Gell-Mann matrices,
again dropping \(h_0\). SU(3) does not close as simply as SU(2): the anticommutator of two generators brings in the totally symmetric tensor
whose conventions follow [MSW68]. Two invariants control the result: the norm \(|h|^2 = h_i h_i\) and the cubic invariant
where \((h \star h)_i = d_{ijk} h_j h_k\) is the star product. These equal \(\mathrm{Tr}(H_0^2)/2\) and \(\mathrm{Tr}(H_0^3)/2\) respectively, with \(H_0\) the traceless part — a relation the test suite checks.
The characteristic equation of \(-H_0\),
has three real roots \(\psi_m\), given in closed trigonometric form. The coefficients of
then follow by Lagrange interpolation over those roots:
Degenerate spectra
The denominator \(3\psi_m^2 - |h|^2\) is the derivative of the characteristic polynomial, so it vanishes exactly at a repeated root. Two degenerate cases are handled separately, and exactly:
\(|h|^2 = 0\), where the Hamiltonian is proportional to the identity and \(U_3 = \mathbb{1}\);
a doubly degenerate root \(\psi_a = \psi_b \neq \psi_c\), where the spectral decomposition collapses onto a single projector,
\[U_3 = e^{i \psi_a L}\mathbb{1} + \left(e^{i \psi_c L} - e^{i \psi_a L}\right) P_c , \qquad P_c = \frac{h_k \lambda^k + \psi_a \mathbb{1}}{\psi_a - \psi_c} ,\]which is linear in \(h_k\) and so needs no matrix algebra.
The argument of the arc cosine that produces the roots lies in \([-1, 1]\) for any Hermitian Hamiltonian, but only up to round-off; it is clipped, so that a marginally out-of-range value cannot yield complex roots and a non-unitary evolution operator.
Four flavors
Everything above generalizes with \(3 \to 4\). Expand in the fifteen generalized Gell-Mann matrices,
again dropping \(h_0\). Three things are new, and each of them is a consequence of SU(4) having rank three where SU(3) has rank two.
A third invariant. The traceless part carries
the first two being the \(|h|^2\) and \(\langle h \rangle\) of the three-flavor case. Taking them from traces means the SU(4) \(d\) tensor — a \(15\times15\times15\) table — is never built.
A quartic, which still solves. The characteristic equation becomes
and Euler’s reduction turns it into the resolvent cubic
whose roots are \(z_i = (\psi_i + \psi_j)^2\) — real and non-negative precisely because \(\tilde{H}\) is Hermitian. So the same trigonometric formula used at three flavors solves it, and then
The SU(3) machinery is literally nested inside the SU(4) solution.
A longer star-product tower. The three-flavor identity \((h \star h) \star h = \tfrac13 |h|^2 h\) is a Cayley-Hamilton accident of \(n = 3\) and is false at \(n = 4\) — over two thousand random Hamiltonians the two sides differ by a median of 56%, with the central 90% of draws between 34% and 116% — so the third rung enters as independent data. The median is the stable statistic here: the extremes of a finite draw are not, and grow with the sample, which is why a range of largest-and-smallest is not quoted.
That also exposes the general-\(n\) pattern: a numerator of degree \(n-2\) in \(\psi_m\), a star tower cut off at length \(n-1\), and always \(\chi'(\psi_m)\) underneath.
Why the method stops at four
Not for want of effort, and not because \(n = 5\) is uninteresting. The whole construction rests on one thing: that the eigenvalues of the traceless Hamiltonian — the roots of its characteristic polynomial — can be written down in radicals, as an explicit formula in the invariants.
That polynomial has degree \(n\). Quadratics, cubics and quartics are solvable in radicals; the Abel-Ruffini theorem says the general quintic is not, and Galois theory says why: the symmetric group \(S_5\) is not soluble, while \(S_2\), \(S_3\) and \(S_4\) are. At \(n = 5\) there is no formula to write, and the shortfall is a theorem rather than a gap in anyone’s algebra.
So the closed-form road ends at four, and it ends for a reason external to neutrino physics entirely.
What does not end there is the philosophy. Nothing above the eigenvalues needs radicals: the interpolation over the roots, the fact that no eigenvectors are ever required, and the whole probability construction go through for any \(n\). Feed numerically computed eigenvalues into the same Sylvester sum and the method degrades gracefully rather than breaking — which is what a general-SU(\(n\)) treatment would do, and what codes carrying SU(\(N\)) expansions to \(n = 6\) in fact do. It would no longer be a closed form, which is this library’s reason to exist, so it is out of scope here.
Four flavors is therefore both the natural stopping point and a useful one: it is exactly what 3+1 sterile scenarios need.
Stiff spectra, and what they cost
First, the perspective, because the numbers below are small enough to be misread as a problem.
None of this is near any measurable effect. Oscillation probabilities are confronted with data at the per-cent level at best, and the systematic uncertainties of a real experiment dominate long before the fourth decimal place. Both figures in this paragraph are errors on a probability. The table further down measures something else — the relative error of the latent roots themselves — so its numbers are smaller and are not comparable to these: an error in a root reaches the probability as a phase error \(\delta\psi\,L\), amplified by the baseline. Even the worst probability here — the unrefined four-flavor result at \(5\times10^{-7}\) — sits four or five orders of magnitude below anything an experiment can resolve, and the refined one at \(10^{-9}\) is far beyond any physics requirement.
So why care? Three reasons, none of them about a single probability:
The claim. This library says it computes probabilities exactly, with no approximation beyond round-off. A probability wrong by \(5\times10^{-7}\) is still round-off-limited in a sense, but it is not the same claim, and the difference should be stated rather than glossed.
Composition.
slabsandearthmultiply evolution operators across many layers, so a per-layer error accumulates. What is invisible in one probability need not stay invisible across a hundred.Regression testing. A suite that pins agreement at \(10^{-9}\) catches a real mistake; one that pins it at \(10^{-6}\) has room for a bug to hide in.
Now the mechanism. A 3+1 Hamiltonian with \(\Delta m^2_{41} \sim 1\) eV2 has a stiff spectrum: the eigenvalues span four orders of magnitude, with three of them clustered. The invariants \(I_2, I_3, I_4\) are sums over that spectrum, so forming them in double precision compresses a \(4\times4\) matrix into three numbers and loses what separates the cluster. Perturbing the three invariants at the \(10^{-16}\) level — their own rounding — moves the roots by \(6\times10^{-11}\) relative at \(\Delta m^2_{41} \sim 1\) eV2, and by more as the splitting grows: the amplification from coefficients to roots measures \(2.3\times10^9\), which on the stiffest Hamiltonian in the reference puts the closed form at \(2.2\times10^{-7}\). That is a property of the problem, not of the solver: it is the classic ill-conditioning of polynomial roots with respect to their coefficients, and it means no better root-finder helps. Deflating the quartic to a cubic first was tried, and does not.
There are two ways out of a bottleneck like that, and the library has shipped both.
Stop asking the invariants. One Newton step on
refines the roots using the Hamiltonian entries at full precision, which
never pass through the three-number bottleneck. A second step changes
nothing — Newton doubles the correct digits, and one already reaches that
route’s floor — so exactly one is taken. This was the whole answer in
1.12.0, and oscprob4nu.POLISH_ROOTS still switches it.
Or widen the bottleneck. The compression is only fatal because
\(10^{-16}\) is what a float64 coefficient carries. Form
\(I_2, I_3, I_4\) in double-double arithmetic — each number a pair
of float64 limbs, together holding some 32 digits — and the same
\(2.3\times10^9\) amplification acts on \(10^{-32}\) instead,
landing at \(10^{-23}\). The roots are then limited by the final
rounding of the answer to float64 rather than by the algebra at all.
One Aberth sweep, also in double-double, carries the quartic’s roots to that
limit. This is the default since 1.13.0, and
oscprob4nu.ROOT_STRATEGY chooses between the two.
Important
Four-flavor results moved in 1.13.0. The roots now come from a
different route, so any probabilities_4nu() can differ
in its last digits from what 1.11.0 returned. Measured on 3+1 spectra at
0.5, 1 and 10 GeV with \(\Delta m^2_{41}\) of 0.1 and 1
eV2, over 1300 km: bit-identical in two of those six, and at
most \(1.2\times10^{-10}\) in the rest.
The new values are the accurate ones — the roots behind them are good to
\(3.6\times10^{-17}\) where the older route reached
\(2.2\times10^{-7}\) — and both sit orders of magnitude below
anything an experiment resolves. Neither setting of
oscprob4nu.ROOT_STRATEGY reproduces 1.11.0 exactly, since
'eigensolver' is LAPACK rather than the closed form 1.11.0 used;
psi_roots_4nu() still solves the quartic in closed form
for anyone who wants that specific computation.
What each gains, and what the alternatives gain
Measured against ground truth from mpmath at fifty decimal digits, worst
of the nine Hamiltonians in tests/stiff_reference.json, which run from
the physical 3+1 range out to \(\Delta m^2_{41} = 1000\) eV2
and down to a pair separated by \(10^{-16}\). Cost is a full
probabilities_4nu() over a 100 000-point stiff stack
through the compiled kernel, relative to what 1.12.0 shipped:
Strategy for the latent roots |
Relative error |
Cost |
Keeps the closed form? |
|---|---|---|---|
Closed form alone |
2.2e-07 |
— |
yes |
Closed form + one Newton step |
3.9e-16 |
1.00x |
yes |
|
6.9e-16 |
0.9-1.0x |
no |
Invariants in double-double |
3.6e-17 |
1.15-1.3x |
yes |
The errors are exact and reproducible to the digits shown; the costs are ranges on purpose. Timed in alternated pairs to cancel machine drift, the double-double row lands at 1.18x, 1.20x and 1.25x by three different methods, while the per-pair spread runs from 0.76x to 1.72x — so a tighter figure would be describing one machine’s load rather than the algorithm. On the same evidence the Newton step in the row above costs under a tenth, which is inside that spread.
The first row has no cost against it because it is no longer a route a
caller can select: psi_roots_4nu() still solves the
quartic in closed form, and its contract is unchanged, but nothing in the
probability path is built on it any more. Its error is quoted on the same
nine Hamiltonians as the rest, which is what it is there to show.
Four things in that table are worth reading twice.
The Newton step is more accurate than LAPACK, by about a factor of two.
That is not a fluke: eigvalsh reduces the matrix by Householder and QR
similarity transforms, each carrying a backward error of order
\(\epsilon \|H\|\), while the Newton step converges onto the root of
\(\det(\psi\mathbb{1} - \tilde{H})\) for the matrix it was handed.
Double-double is the only row that is not fighting the conditioning.
Every other entry accepts a \(10^{-16}\) coefficient and then tries to
recover from it, whether by refining against the matrix or by avoiding the
polynomial altogether; this one denies the amplification anything to work
on. It is also the only row whose accuracy is set by the output format
rather than by the method, which is why more Aberth sweeps do not improve
it: one, two and three give bitwise identical answers on all nine
Hamiltonians. That one suffices is a property of the start, not of the
iteration — the eigensolver’s quartet with the residual-trace shift
removed in double-double rather than in float64. Removing it in
float64 throws away the low limb, puts the start an ulp out, and costs
a second sweep to recover; the count was two for exactly as long as that
was how it was done.
Extended precision is still a poor trade, and instructively so. Solving
the same quartic in numpy.longdouble was measured at 4.5e-11 — under
one extra digit, rather than the three its wider mantissa suggests, because
the cluster amplifies whatever coefficient error survives — while being
slower for not being hardware-vectorised, and silently float64 on Apple
Silicon and on Windows, where the “fix” would quietly do nothing at all.
Double-double gets four orders more than that, portably, out of pairs of the
float64 every machine has.
The double-double route does call eigvalsh, once per element, and only
for a starting quartet. That is a smaller concession than it looks:
what a start has to supply is not accuracy but four distinct basins, and its
own error is discarded by the first sweep. Euler’s closed form is twice as
fast and exact on every stiff case in the reference, and still cannot be
used, because on a cluster Aberth converges linearly at ratio one half —
from \(10^{-7}\) five sweeps reach \(3.8\times10^{-9}\), and it
would need some thirty. A backward-stable Hermitian eigensolver separates a
cluster as well as float64 allows, in one call. Durand-Kerner was
measured in place of Aberth as well, at one division per root against five,
and rejected for being non-monotone in the sweep count: 3.9e-16, then
9.7e-17, then 1.9e-16.
The honest summary on cost is unchanged in shape: four flavors costs more per point than three, and the accurate routes cost more than the fast ones.
Why the refinement is not applied selectively
The obvious saving is to refine only the elements that need it, the way
oscprob3nu.SMALL_BATCH and fastkernels.MIN_BATCH dispatch on
a measured threshold. It was measured, and it does not work. The result is
recorded here so that it is not rediscovered.
Two criteria were tried, on 6300 Hamiltonians spanning clustered, doubly
paired and generic spectra, against eigvalsh. A criterion is safe at a
given cut only if every element below the cut is more accurate than the
target; the question is how many elements a safe cut can skip.
The gap-based amplification \(\max_m |\psi|^3_{\max}/|\chi'(\psi_m)|\), which is what perturbation theory suggests, since the root sensitivity goes as \(1/\chi'\). It predicts the error well for a single cluster and badly for two degenerate pairs, a family it does not model: those reach \(1.7\times10^{-10}\) at an amplification of ten, where the criterion expects round-off. The largest safe cut is about 2.3, which is close to the smallest value the indicator ever takes — so it skips nothing.
A matrix residual, comparing \(\prod_m \psi_m\) against \(\det \tilde{H}\), which is \(\chi\) evaluated at zero and costs one determinant instead of four. Being built from the matrix it cannot lose information the way a gap heuristic does, but it is one scalar constraint on four roots, and errors cancel in the product: there are samples with a residual of \(10^{-17}\) and a root error of \(5\times10^{-5}\). It also skips nothing safely.
The second failure points at the general reason. A criterion complete enough to certify all four roots has to evaluate \(\chi\) at all four roots — and that is the refinement. The check and the fix are the same computation, so there is nothing to save by doing the check first.
That conclusion held up when it was retested for the double-double route, where it cuts the other way and hands over a free certificate. An Aberth sweep evaluates \(\chi\) at four roots by construction, so the size of the step it takes says whether the answer has converged — around \(10^{-32}\) relative once it has, against \(10^{-9}\) while still crawling up a cluster, a gap wide enough that any threshold between them works. That made a genuinely cheap arrangement possible: start from Euler’s closed form, which costs nothing, and fall back to the eigensolver only for the elements the certificate rejects. It was built and measured, and dropped — the fallback fires on 80% of stiff spectra and 99% of those with a tight pair, which leaves it slower than simply starting from the eigensolver every time, for twice the code. The criterion was sound; the saving it was guarding was not there.
Note also who would benefit. The four-flavor module exists mainly for 3+1, and a 3+1 scan is stiff at every point, so even a working criterion would skip nothing on the workload that motivates the module, while adding its own cost. Unconditional refinement is therefore not a compromise: it is what the measurement supports.
Finally, this is specific to four flavors rather than a general caveat. The
same measurement on oscprob3nu gives \(10^{-14}\), because there
\(\Delta m^2_{31}/\Delta m^2_{21}\) is 34 rather than 13500.
That figure has since been reached independently. Compared against
nuSQuIDS, which integrates the
density matrix numerically and shares nothing with this expansion, the
four-flavor probabilities agree to \(4\times10^{-16}\) where the spectrum
is benign and \(3\times10^{-10}\) where it is stiffest — the latter
being this limit, arrived at from outside rather than asserted from within.
Notebook 17
attributes that residual: against scipy.linalg.expm on the same
Hamiltonian our error is the same size, so it is ours and not a
disagreement.
Degenerate spectra at four flavors
Reconstructing \(U_4\) from its roots divides by their differences, so a repeated root needs care. Rather than branch on a tolerance, the exponential is interpolated over the roots in Newton form, with divided differences: a repeated node is then a derivative, and for the exponential that derivative is known exactly, \(f^{(k)}(\psi)/k! = (-iL)^k e^{-i\psi L}/k!\).
The alternative — solving the Vandermonde system for the Cayley-Hamilton coefficients — is singular the moment two roots coincide, which includes a Hamiltonian proportional to the identity and any triply degenerate spectrum. The Newton form handles every degenerate case with no special branch at all.
Sign conventions
The vacuum Hamiltonians are built so that a positive matter potential added to the \(ee\) entry describes neutrinos:
with \(M^2 = \mathrm{diag}(0, \Delta m^2_{21}, \Delta m^2_{31})\) and \(U\) the PMNS matrix in the standard PDG parametrization.
This is worth stating explicitly because an overall sign flip of the vacuum
term alone is invisible in vacuum — for a real Hamiltonian the
probabilities are invariant under \(H \to -H\) — and yet reverses the
sign of the matter potential relative to it, turning neutrinos into
antineutrinos and moving the Mikheyev-Smirnov-Wolfenstein resonance to the
other side. For antineutrinos, pass a negative VCC (and conjugate the
CP phase).
Ordering of the probabilities
oscprob3nu.probabilities_3nu() returns the nine probabilities with the
initial flavor varying slowest, so the returned tuple reads
with \(P_{\alpha\beta} \equiv P(\nu_\alpha \to \nu_\beta) = |[U_3]_{\beta\alpha}|^2\). Note the index order: the evolution operator is indexed (final, initial), the probabilities (initial, final).
Cost
The \(d\) tensor is constant and is tabulated once at import time as a dense \(8\times8\times8\) array; the star product and the two invariants are contractions against that table.
A single three-flavor probability evaluation takes about eight microseconds, and a two-flavor one about one. For scans, pass arrays rather than looping: the routines accept a stack of Hamiltonians, a stack of baselines, or both, and evaluate the stack in one pass. Measured against the equivalent Python loop, on 2000 points:
Scan |
Speedup |
Also comparable to |
|---|---|---|
Versus baseline (one \(H\), many \(L\)) |
~21x |
one |
Versus energy (many \(H\), one \(L\)) |
~23x |
batched |
Oscillogram, 100 x 100 |
~37x |
|
Two flavors, versus baseline |
~99x |
These are the ratios of the same four measurements tabulated on the
landing page and in README.md, which state them as
absolute timings; the numbers here are those timings divided. They used to
be quoted independently and had drifted apart — 30x against 21x, and 70x
against 99x — which is why they are now derived from one set of figures
and guarded together by tests/test_documented_figures.py.
These ratios have narrowed across successive releases even as both sides got quicker: the scalar path has itself sped up several-fold, so the loop being compared against is no longer as slow as it was. In absolute terms the vectorised scan is faster than it has ever been.
The two scans differ because the latent roots depend on the Hamiltonian alone. Scanning one Hamiltonian over many baselines solves the characteristic equation once and then only evaluates \(e^{i\psi_m L}\), so almost all the work is amortised; scanning over energy changes the Hamiltonian at every point and must solve it each time.
The vectorised path is as fast as diagonalising with LAPACK, which is the
honest comparison to draw: the SU(3) route’s advantage is that it is a
closed form, not that it outruns eigh. What the vectorisation removes
is the disadvantage it used to carry.
Short stacks
A batched call carries a couple of hundred microseconds of fixed cost
whatever its length — allocating and reducing a dozen small arrays — so
for a handful of points it spends more on the machinery than the scalar path
spends on the whole job. Stacks below
oscprob3nu.SMALL_BATCH are therefore evaluated one element at a time.
The thresholds are measured, not guessed: thirteen elements for three flavors, twelve for two. They sit close together even though the two-flavor expansion does much less work per element, because what has to be amortised is the array machinery’s fixed cost rather than the arithmetic. Nothing about this is visible from the outside; the answers are the same either way.
The compiled backend
With fastkernels — that is, with Numba, which is installed with the
package as of 1.13.0 — the batched paths are compiled instead. The NumPy path evaluates the expansion as a
chain of whole-array operations, so a stack of \(N\) Hamiltonians makes
roughly fifteen passes over \(N\)-element arrays, each writing a
temporary that the next reads back. The compiled kernel does the same
arithmetic one element at a time, keeping every intermediate in registers,
and spreads the elements over the available cores. Against the NumPy path:
Stack |
Speedup |
|---|---|
200 000 energies, four flavors |
~19x |
20 000 energies, four flavors |
~18x |
200 000 energies, three flavors |
~15x |
20 000 energies, three flavors |
~9x |
100 x 100 oscillogram |
~3.5x |
200 000 baselines, two flavors |
~1.5x |
Four flavors gains the most, and the reason is worth stating, because it is not that the kernel is cleverer there. That expansion needs a quartic, a Newton refinement of its four roots against the matrix, and a Newton-form reconstruction; as whole-array operations that is some forty passes rather than fifteen, one of them a batched \(4\times4\) determinant. Done one element at a time none of it leaves the registers, so the path carrying the most fixed cost is the one with the most to shed.
The two-flavor row is the honest caveat: that path 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 — for a scan over baselines, the same matrix repeated, which costs 2.5 ms to copy at two hundred thousand points.
So the backend is not used unconditionally. fastkernels.worthwhile()
declines it below the per-flavor thresholds in
fastkernels.MIN_BATCH, which were found by alternating the two paths
and taking the best of nine rounds each: for three and four flavors the
kernel wins at every size, for two it wins from about fifty thousand
elements. A backend
that is sometimes slower than the path it replaces is worse than no backend,
and a test asserts the thresholds are honoured.
The scalar path is deliberately left uncompiled. One probability takes about eight microseconds; compiling it would save most of that, at the cost of a multi-second pause on a user’s first call.
Checking the input, and what it costs
The expansion assumes a Hermitian Hamiltonian, and one that is not Hermitian does not fail loudly: the probabilities it returns still sum to one, so the check a caller would actually apply cannot tell that the answer is meaningless. Every entry point therefore verifies it, and the tolerance is relative to the largest entry, so a matrix assembled in floating point passes — everything the sample-Hamiltonian modules build is Hermitian to about \(2 \times 10^{-17}\) relative, against a tolerance of \(10^{-12}\).
The cost is stated rather than buried, because it is larger than it looks. Validating a stack is a pass over it, which is the same order of work as evaluating it, and the compiled kernel has made evaluating it fast; on a large stack the check therefore dominates. Interleaved, best of fifteen rounds each:
Stack |
Two flavors |
Four flavors |
|---|---|---|
2 000 points |
1.5x |
1.3x |
200 000 points |
5.7x |
3.2x |
Three flavors sits between them, at 1.8x and 3.9x. Two ways of making it
cheaper were measured. Comparing real and imaginary parts separately, rather
than forming \(H - H^\dagger\), avoids a temporary the size of the stack
and a square root per element, and is what the code does. Replacing
np.abs(...).max() with reductions that allocate nothing was the obvious
next step and came out 1.4x slower, because those views are strided over
the complex array while np.abs reads it contiguously.
It defaults to on regardless, because a library that silently returns
meaningless numbers costs its user more than the check does.
oscprob3nu.CHECK_HERMITICITY turns it off, per module.
Degenerate spectra on the batched path
The degenerate branch cannot be taken elementwise inside a vectorised expression. The general formula is therefore evaluated everywhere, with vanishing denominators replaced by one, and the affected elements are then recomputed individually. Degeneracy is measure-zero among floating-point Hamiltonians, so that fallback loop is empty in essentially every real use.