magnus

The Magnus paper is out!

Standard

I have posted the paper that accompanies Magνs, the code I released in August for computing neutrino oscillation probabilities.

In vacuum and in matter of constant density, the oscillation probability has an exact formula. Where the density changes along the path, as through the Earth, the Sun, or a supernova, it does not, and the equation of motion has to be integrated. The usual approaches both become slow at high accuracy. Cutting the matter profile into a staircase of constant-density steps gives an error that falls only as the square of the step width. Handing the equation to a general-purpose solver means resolving every oscillation, so the cost grows with the number of oscillations along the path, which can reach millions for a neutrino leaving the Sun. On top of that, each new experiment or model often means writing and validating a new solver.

Magνs uses the Magnus expansion instead. Over a stretch of the path within which the density may vary, it writes the evolution of the neutrino as a single exponential built from integrals of the Hamiltonian. The error falls as up to the tenth power of the stretch’s length, and the cost is set by how fast the density changes, not by how many times the neutrino oscillates. Three advantages follow, whatever the Hamiltonian contains.

Robust. The evolution is exactly unitary at any order and any accuracy setting, so the probabilities always add up to one, to round-off.

Fast. A typical probability takes a few milliseconds, and a scan over energy or arrival direction is a single batched call, one to two orders of magnitude cheaper per point than computing the points one at a time.

Flexible. The Hamiltonian is any Hermitian matrix, at any number of flavors. Non-standard interactions, sterile neutrinos, Lorentz-invariance violation, pseudo-Dirac neutrinos, or an interaction nobody has written down yet all go through the same call.

Magνs also chooses its own method. You state a tolerance, and it decides how finely to divide the path. Where a cheaper route exists, it takes it, such as following the slowly changing eigenstates of the Hamiltonian across the Sun. Besides oscillating probabilities, it returns phase-averaged ones, the quantity solar and astrophysical neutrino experiments measure. The aim is to move the effort of studying a new experiment or model from writing a solver to specifying its physics.

The paper works through examples across the field: the Earth’s layers and core, the Sun, a supernova shock front, geoneutrinos, the flavor composition of astrophysical neutrinos, a jet inside a collapsing star, turbulent matter, and a Hamiltonian supplied by the user. Each is a variation on one call.

Install Magνs via PyPI:

pip install magnuspy

Documentation: https://mbustama.github.io/Magnus/

Notebooks: https://mbustama.github.io/Magnus/tutorials.html

Read more at:

Magνs: neutrino oscillation probabilities for any Hermitian Hamiltonian, any number of flavors, and any matter profile
Mauricio Bustamante
2610.07159 hep-ph

Magnus: neutrino oscillation probabilities for any Hamiltonian, in any number of flavors

Standard

Documentation: https://mbustama.github.io/Magnus/

Notebooks: https://mbustama.github.io/Magnus/tutorials.html

Source code: https://github.com/mbustama/Magnus

PyPI: https://pypi.org/project/magnuspy/

Oscillations in matter are usually computed by pretending the matter holds still: chop the path into slabs of constant density, solve each one exactly, multiply. It is a good trick, and for a neutrino crossing the Earth it is nearly the truth. For the Sun, for a supernova, for anything with a front moving through it, the pretending is where the physics is.

Magνs is my new package for that case. It computes oscillation probabilities for an arbitrary number of flavors, under any Hamiltonian, whether or not it depends on time. Rather than stepping the Schrödinger equation forward, it integrates the Hamiltonian across each slab and exponentiates the result — the Magnus expansion. The useful consequence is structural: however hard the expansion is truncated, the truncation stays inside the Lie algebra, so the evolution operator is exactly unitary by construction. Probabilities come out non-negative and summing to one at machine precision, at every accuracy setting, not only the expensive ones.

What it makes easy

  • A density that actually changes — a tabulated solar model, a supernova profile from a simulation, a shock front.
  • An accuracy instead of a slab count — ask for a tolerance and the refinement finds the resolution for you.
  • Any number of flavors — two through five, sterile states (3+1, 3+2) included.
  • The phase-averaged probability, in closed form — the quantity a solar experiment actually measures, returned directly rather than reconstructed by averaging a scan.
  • Neutrinos through the Earth — you give a direction, not a profile.

Install with:

pip install magnuspy

The import package is magnus. Python 3.10 or newer, GPL-3.0-only, and 27 example notebooks that are executed in CI, so the numbers printed in them cannot quietly drift.

A minimal example

A 10 GeV muon neutrino arriving from below, crossing the Earth:

import magnus.oscprob as oscprob
import magnus.globaldefs as gd
import magnus.earth as earth
# A 10 GeV muon neutrino crossing the Earth, arriving from below.
costhz = -0.4
L = earth.distance_traveled_inside_earth(costhz)*gd.UNIT_KM   # chord: km -> eV^-1
P = oscprob.osc_prob_3nu_earth(10.0*gd.UNIT_GEV, costhz=costhz, L=L)
print('Pme = %.5f, Pmm = %.5f, Pmt = %.5f'
      % (P[gd.NUMU][gd.NUE], P[gd.NUMU][gd.NUMU], P[gd.NUMU][gd.NUTAU]))

which prints

Pme = 0.10273, Pmm = 0.01165, Pmt = 0.88562

The oscillation parameters are never mentioned, because they default to the current global fit (NuFIT 6.1). Neither is the Earth’s density: the wrapper walks the PREM profile for the chord that the zenith angle picks out. What you supply is a direction and an energy.

How it differs from NuOscProbExact

Both packages are mine, and they are not competitors — the honest way to put it is as a boundary rather than a winner. NuOscProbExact solves each constant-density slab in closed form, through SU(2), SU(3) and SU(4) expansions. It is exact up to round-off, there is no slab count to argue about, and it is fast. Where a closed form exists and the accumulated phase is large, use the closed form — including the Earth through PREM at three flavors, where NuOscProbExact is about 20× cheaper per call than Magνs. That is its home ground.

Magνs integrates across each slab instead, which is exactly what lets it follow a density that changes while the neutrino is still inside it. Reach for it when the profile varies continuously, when you would rather ask for an accuracy than guess a resolution, when you need five flavors — the SU(N) closed forms stop at SU(4) — or when the Hamiltonian is something nobody has diagonalized.

Two measurements make the difference concrete. On a smooth profile, multiplying constant-density slabs together bottoms out near 2.5 × 10⁻¹¹ and then gets worse: past roughly 16 000 slabs, round-off from composing that many matrix products outgrows anything a finer grid buys. Magνs keeps going to 2.9 × 10⁻¹³. And for the Sun, it returns 40 phase-averaged energies in about 0.7 s, where nuSQuIDS needs on the order of ten minutes merely to reach the solver tolerance at which its output is a probability at all.

The two packages deliberately share conventions, units and parameter defaults, so a calculation can be moved from one to the other and used as a cross-check on itself. I do this often.

Some physics it’s meant for

Each of these has a worked notebook in the tutorial gallery, and there is a systematic comparison against other codes in the documentation.

  • Long-baseline beams and atmospheric neutrinos through the Earth.
  • Solar neutrinos on a real tabulated model, including the averaged observable.
  • Supernova neutrinos across a shock front — where the width of the front is what decides which method the problem belongs to.
  • Sterile states (3+1 and 3+2), non-standard interactions, Lorentz-invariance violation.
  • A Hamiltonian of your own — the escape hatch, and the reason the rest of the list is not exhaustive.