Phase-averaged probabilities

This page documents four keywords and the module that implements them, magnus.avgprob. The keyword average asks for the phase-averaged probability; average_spread sets the width of the average; average_initial_state sets the state the neutrino starts in (see The initial state); and average_n_samples sets the number of energies sampled where there is no closed form (see When there is no closed form). Every scenario function and every osc_prob_{2,3,4,5}nu_* wrapper accepts them. For a varying profile, the average reuses the machinery of Adiabatic + Magnus hybrid strategy.

Tip

In short. average=True returns the probability averaged over an energy spread, 10% by default (average_spread). When every oscillation phase is large — the Sun, astrophysical sources — the result does not depend on the spread. When some phase is not — a beam experiment at its first oscillation maximum, for instance — it does, and a PhaseAveragingWarning says so: set average_spread to the energy resolution of the measurement. average_spread is the width of a Gaussian in \(\ln E\), except on a profile with declared discontinuities, where it is the half-width of a top-hat energy window (When there is no closed form).

The problem: a phase nobody can resolve

A neutrino from an astrophysical source arrives with an oscillation phase

\[\Delta \phi = \frac{\Delta m^2 L}{2E} \simeq 2.53 \times \frac{\Delta m^2/\text{eV}^2 \times L/\text{km}} {E/\text{GeV}}\]

of order \(10^{16}\) for a TeV neutrino from 100 Mpc away. No ingredient of that number is known to anything close to the precision the phase would demand: not the distance, not the size of the source and not the energy of the neutrino. Whatever the true phase is, the measurement integrates over many complete cycles of it.

Computing such a probability by propagation is doubly unattractive. It is expensive — resolving \(10^{16}\) radians is exactly the regime that defeats slab refinement — and it is pointless, because every oscillatory term is averaged away by the integration that the measurement performs.

The averaged limit

Write the amplitude in the basis that diagonalizes the Hamiltonian, \(\mathbb{H} = \mathbb{V} \,\mathrm{diag}(\lambda_i)\, \mathbb{V}^\dagger\):

\[A(\nu_\alpha \to \nu_\beta) = \sum_i \mathbb{V}^*_{\alpha i} \mathbb{V}_{\beta i}\, e^{-i \lambda_i L} .\]

The probability \(|A|^2\) contains a diagonal part and interference terms carrying \(e^{-i(\lambda_i - \lambda_j)L}\). Averaging over the phase leaves only the terms whose phase does not vary:

\[\boxed{\;P(\nu_\alpha \to \nu_\beta) = \sum_i |\mathbb{V}_{\alpha i}|^2\, |\mathbb{V}_{\beta i}|^2\;}\]

This is the exact \(L/E \to \infty\) limit, reached when every relative phase runs through many cycles across what the measurement cannot resolve, and it costs one matrix product rather than an integration. The limit has three properties, none of which holds away from it:

The result is symmetric in \(\alpha \leftrightarrow \beta\), so the averaged probability is the same in both directions.
In vacuum, it is identical for neutrinos and antineutrinos, since the antineutrino replaces \(\mathbb{V}\) by \(\mathbb{V}^*\) and \(|\mathbb{V}^*|^2 = |\mathbb{V}|^2\). The phase \(\delta_{\rm CP}\) still enters through the magnitudes \(|\mathbb{V}_{\alpha i}|\), but produces no difference between the two. In matter, they differ because the potential changes sign. A Lorentz-violating term of even n_liv also changes sign, so the two differ in vacuum as well.
It does not depend on the baseline, and depends on the energy only through the eigenvectors. In standard vacuum oscillations, those do not change with energy, since \(\mathbb{H} \propto 1/E\), so the result is a constant and a single matrix serves an entire flux calculation. In matter, it changes with the energy, through the competition between the potential and the vacuum term.

The phase average

A measurement averages over what it cannot resolve, and the limit assumes that this includes every relative phase. It need not. A detector with an energy resolution of 10% averages a phase over the range the phase covers across that resolution. A phase that barely changes with energy is not averaged at all. Over 1000 km at 1 GeV, the atmospheric phase is about 6 rad, and the limit misses the oscillation entirely.

Instead, average=True returns the phase average. Every interference term keeps its phase at the central energy and is weighted by how much that phase varies across a relative energy spread \(\sigma\):

\[P(\nu_\alpha \to \nu_\beta) = \sum_{ij} \mathbb{V}^*_{\alpha i} \mathbb{V}_{\beta i} \mathbb{V}_{\alpha j} \mathbb{V}^*_{\beta j}\, e^{-i\phi_{ij}}\, e^{-\sigma^2 \phi_{ij}'^2/2} , \qquad \phi_{ij} = (\lambda_i - \lambda_j) L , \quad \phi'_{ij} = \frac{d\phi_{ij}}{d\ln E} .\]

The factor \(e^{-\sigma^2\phi'^2/2}\) is the average of \(e^{-i\phi}\) over a Gaussian spread of width \(\sigma\) in \(\ln E\), with the phase taken to first order in the offset. In vacuum, \(\phi' = -\phi\), and at \(\sigma = 10\%\) the fraction of an interference term that survives is

phase [rad]	1	2π	10	20	30	40	50
surviving	0.995	0.82	0.61	0.14	0.011	3e-4	4e-6

Hence, the phase average is the boxed limit where every phase runs through many cycles, the oscillation probability itself where none does and a smooth weighting in between. The spread is average_spread, 0.1 by default (magnus.avgprob.AVG_PHASE_SPREAD), and it should be the resolution of the measurement.

Only the phases are averaged. Mixing, the crossing amplitudes of the next sections, and the eigenbases at the two ends of the path stay at the central energy, so a probability without interference is returned unchanged. Averaging the whole probability over energy would smooth those too: on the 'BS05-AGS-OP' solar curve from 0.1 to 20 MeV it moves 26 of 40 energies by more than 1e-4, and by 5.2e-4 at most, whereas the phase average moves none.

The slopes come from the Hellmann–Feynman theorem, \(d\lambda_i/d\ln E = \langle v_i|\, d\mathbb{H}/d\ln E\, |v_i\rangle\), with \(d\mathbb{H}/d\ln E\) a central difference of the Hamiltonian in \(\ln E\) with step 1e-3: two more evaluations of \(\mathbb{H}\) per energy, and no further eigendecomposition. A slope smaller than the round-off of that difference is replaced by minus the phase of its pair, its value in vacuum. Without that, a pseudo-Dirac pair split by 1e-21 eV² reads a slope of 0.2 rad at 100 TeV over 100 Mpc, against a true 8e-5 rad, and its result depends on the spread although the pair is coherent.

A Hamiltonian that does not depend on energy — a matrix, or a function of position alone, passed to osc_prob_energy_baseline (see Usage) — has no slope for a spread to act on. average=True then averages the pairs of levels whose phase has grown large and keeps the others coherent: at a short baseline, where no phase has grown, it returns the coherent probability. Away from the limit, the three properties above fail: \(P_{\alpha\beta} \neq P_{\beta\alpha}\) in general, and CP violation survives in the interference terms that remain.

Every point is first computed in the limit. The limit is returned, bit for bit, wherever the phase average agrees with it to within 1e-4, or to within the tighter of rtol and atol if smaller, both absolutely and relative to each probability.

Coherence, and where the spread matters

Whether a pair of eigenvalues has decohered is a statement about that pair, not about the spectrum as a whole. The phase average decides it per pair through the weight; the functions of magnus.avgprob that return the limit decide it with a threshold. magnus.avgprob.averaged_probabilities_constant_hamiltonian() groups the spectrum into blocks of mutually coherent eigenvalues (pairs whose phase is below magnus.avgprob.DECOHERENCE_PHASE_THRESHOLD, \(2\pi\)), and sums coherently inside each block,

\[P(\nu_\alpha \to \nu_\beta) = \sum_{b} \Big| \sum_{i \in b} \mathbb{V}^*_{\alpha i} \mathbb{V}_{\beta i} \Big|^2 ,\]

which reduces to the boxed expression when every block is a singleton. The block form matters: for a sterile state with a small \(\Delta m^2_{41}\), or any exactly degenerate spectrum, the naive sum is wrong without any warning. With all eigenvalues equal, nothing oscillates and the correct answer is the identity, but the naive sum returns a mixture.

Pseudo-Dirac neutrinos are the case the block form was written for. Each mass eigenstate that carries a sterile partner splits into two states, separated by a \(\delta m^2\) many orders of magnitude below \(\Delta m^2_{21}\). Over an astrophysical baseline, the standard phases have averaged away while every pair is still coherent. Each pair is then one block, and the naive sum is wrong by a factor of two. The phase average needs no blocks for them: a pair whose phase barely moves with energy keeps a weight near one. magnus.hamiltonians.hamiltonians_pseudodirac builds those Hamiltonians, with the pairing selectable per mass state; notebook 29 increases the splitting through the three regimes below, using magnus.avgprob.coherence_report().

For them, use magnus.oscprob.osc_prob_pseudo_dirac_vacuum() with average=True rather than passing the Hamiltonian to osc_prob_energy_baseline(). The generic route takes each pair’s phase from the difference of two eigenvalues, which an eigendecomposition resolves only to about \(10^{-16}\) of the largest, so a splitting far below the standard ones is partly lost. Against a 50-digit reference, over pair phases of 0.5 to 30 rad, it is off by \(2\times10^{-7}\) at \(\delta m^2/\max m^2 = 10^{-10}\), \(7\times10^{-5}\) at \(10^{-12}\), \(6\times10^{-3}\) at \(10^{-14}\) and 0.25 at \(4\times10^{-16}\). The pseudo-Dirac function builds the eigensystem from the pairing instead (magnus.avgprob.phase_averaged_probabilities_pseudo_dirac()): the phase between the two members of a pair is \(\delta m^2_j L/2E\), exactly, and its slope in \(\ln E\) is \(-\phi\). It is within \(1.1\times10^{-16}\) of the reference at every one of those ratios. The measurement is in docs/dev/measurements/issue165_pseudo_dirac_resolution/.

The same per-pair phase decides whether the limit applies at all. A pair is in one of three regimes:

The coherent, intermediate and decohered regimes of a pair of eigenvalues

far below magnus.avgprob.COHERENCE_PHASE_THRESHOLD, the relative phase has barely advanced: the pair is coherent and there is nothing to average;
far above magnus.avgprob.DECOHERENCE_PHASE_THRESHOLD, the cross term has averaged away and the boxed expression is exact;
in between, neither statement holds, and the limit does not describe the result.

magnus.avgprob.coherence_report() names the pairs in that middle band. For the phase average, it is the band where the answer depends on \(\sigma\), and that is what magnus.oscprob.PhaseAveragingWarning reports: it fires where \(|\sigma\, \partial P/\partial\sigma|\) exceeds magnus.avgprob.PHASE_SPREAD_SENSITIVITY_THRESHOLD, 1e-3, the default tolerance. The returned probability is then the average over the requested spread. Asking for the average at a 1000 km baseline does exactly this: at 1 GeV, \(P_{\mu\mu}\) is 0.91 at a 10% spread and 0.97 at 5%.

The initial state

The phase average depends on the state the neutrino starts in, which average_initial_state sets:

'flavor', the default: the flavor state \(\nu_\alpha\), a coherent superposition of the eigenstates at the start of the path, as for a neutrino produced in the medium — a beam, an atmospheric neutrino, a neutrino made in the Earth or in the Sun.
'decohered': an incoherent mixture of those eigenstates, with weights \(|\mathbb{V}_{\alpha i}(l_0)|^2\), as for a neutrino that lost its coherence before reaching the start of the path — one from a distant source crossing the Sun, as in the paper’s solar tomography.

The two differ only by the interference between eigenstates present at the start, each term damped by the spread of its phase from there on, so they agree wherever those phases are large. On a constant Hamiltonian, the flavor start gives the phase average above and the decohered start gives the limit exactly, whatever the phases. On a smooth profile with no non-adiabatic window, a flavor start has a closed form. Before computing it, magnus.avgprob estimates the phase slopes at 17 points; where every interference term would be damped regardless, it returns the decohered result, bit for bit. On a profile with declared discontinuities, the average is over an energy window, starting in a flavor state, and 'decohered' raises ValueError.

import numpy as np
import magnus.oscprob as oscprob
import magnus.globaldefs as gd

osc = gd.load_nufit_params('NuFIT 6.1')
E, L = np.array([1.0*gd.UNIT_GEV]), 300.0*gd.UNIT_KM
for start in ('flavor', 'decohered'):
    P = oscprob.osc_prob_3nu_matter_constant_density(E, L=L, rho=3.0,
        density_matter_is_in_g_per_cm3=True, nu_i=gd.NUMU, nu_f=gd.NUMU, average=True,
        average_initial_state=start, **osc)
    print(start, float(np.ravel(P)[0]))      # 0.3626 and 0.4532

Position-dependent Hamiltonians

When the Hamiltonian varies along the trajectory, there is no single eigenbasis to decohere in. Where every phase along the path has averaged away, the neutrino is carried along the levels of the instantaneous Hamiltonian from the eigenbasis at \(l_0\) to the eigenbasis at \(l_1\), and whatever state it started in, the result is

\[P(\nu_\alpha \to \nu_\beta) = \sum_{ij} |\mathbb{V}_{\alpha i}(l_0)|^2 \, P^{\rm cross}_{ij} \, |\mathbb{V}_{\beta j}(l_1)|^2 ,\]

the standard MSW-plus-decoherence result, generalized here to any number of levels and any number of crossings. The matrix \(P^{\rm cross}\) holds the probability of ending on level \(j\) having started on level \(i\). Adiabatic evolution keeps a neutrino on its level, so \(P^{\rm cross}\) is the identity wherever the adiabatic approximation holds, and departs from it only across a non-adiabatic window.

Those windows are located with the Hellmann–Feynman diagnostic of magnus.adiabatic (see Adiabatic + Magnus hybrid strategy), and the transfer across each one is computed with that module’s own convergence-checked Magnus patch — not with a Landau–Zener formula. The result is an exact treatment of the crossing rather than an asymptotic approximation to it. As a check, the computed hop probability reproduces the analytic Landau–Zener value \(\exp(-2\pi\epsilon^2/|d\Delta/dl|)\) to a few parts in a thousand for a linear crossing, with nothing in the implementation assuming that formula:

Coupling \(\epsilon\)	Computed hop probability	Landau–Zener	Difference
\(10^{-4}\)	0.9382	0.9391	0.1%
\(3\times10^{-4}\)	0.5665	0.5680	0.3%

The search for windows runs on a grid of 200 probes. A feature narrower than their spacing falls between two probes, and no window opens there. To close that gap, Magνs first checks the profile for features that are this sharp and that, if crossed instantly, would move more than magnus.avgprob.SUDDEN_TRANSFER_THRESHOLD (the default tolerance) of probability between levels. Where there is one, the windows are taken from magnus.adiabatic.hybrid_propagator(), which refines its search until it certifies. Of 24 fronts on a supernova shock ray, 0.07–2000 km wide, 16 then agree to within 0.01 with a decohered reference, and the other 8 warn. Where the feature is a discontinuity that no refinement resolves, the call issues magnus.oscprob.UnmarkedDiscontinuityWarning; declare it with t_breakpoints, which selects the energy-window average described below. On every solar profile measured, tabulated models included, the check finds nothing and does not alter the result.

The phase average on a profile

The expression above composes the crossings as probabilities and reads the result out in the eigenbasis at \(l_1\), so it keeps no interference at all. Between two crossings, or between the last crossing and the end of the path, every phase is dropped, however short it is. On a chord through the edge of the Sun at 10 TeV, the last window ends where the path does, the phase that follows is zero and the expression gives 0.40, whereas the true probability at exit is 0.54. magnus.avgprob.averaged_probabilities_adiabatic(), which returns it, reports the pairs for which that composition is unsafe.

The same definition of the phase average extends to a varying profile. An energy offset \(u = \delta\ln E\) moves every instantaneous eigenvalue by \(u\, d\lambda_i/d\ln E\) and leaves the eigenvectors alone. The result is the Gaussian average over \(u\) of the evolution under \(\mathbb{H} + u\, D_{\rm diag}\), with \(D_{\rm diag}\) the part of \(d\mathbb{H}/d\ln E\) diagonal in the instantaneous eigenbasis. magnus.avgprob.phase_averaged_probabilities_adiabatic() computes it without sampling energies along the path. Each non-adiabatic window is an amplitude matrix, computed with the same Magnus patch at a few Gauss–Hermite nodes in \(u\). How fast the phases inside the window change with energy sets the number of nodes, at most 31; beyond that, the nodes form a uniform grid. Each adiabatic stretch between windows is a diagonal phase with an exact slope in \(u\), carried with the parallel-transport phase of its eigenvectors. The density matrix is carried as terms labeled by accumulated slope, whose Gaussian average is analytic, and a term is dropped only once no later stretch can bring its slope back.

Because the definition reaches inside the windows, the answer does not depend on where they are drawn: one window over a stretch, and two windows with the stretch between them, return the same number. Where there is no window, a decohered start carried adiabatically has no interference to keep, and average=True returns the expression above, bit for bit. A start in flavor, the default, keeps the interference between the eigenstates present at production in closed form (see The initial state). For a neutrino produced in the solar core, that interference is damped away entirely, and every solar MSW curve is again the expression above, bit for bit. Transfer between levels outside the windows is neglected, as in any adiabatic calculation. The limit drops the interference such a transfer carries and the phase average keeps it, so the windows are searched at an adiabaticity threshold of 0.01 (magnus.avgprob.PHASE_AVERAGE_WINDOW_THRESHOLD) rather than 0.1. Against a brute-force average of the same definition on five solar chords from 10 GeV to 10 TeV, the result is within 4.2e-5; on the same chords, the limit is off by up to 0.14.

The call’s rtol and atol set the tolerance: the window patches and the stretch phases converge to the tighter of the two, 1e-3 by default. The limit is returned only where it agrees with the phase average, by the rule in the phase-average section. On twenty chords through the solar core, from 30 GeV to 3 TeV, the default moves the probability by at most 4.6e-6 from its value at 1e-5, at 2.4 times the speed. Called directly, magnus.avgprob.phase_averaged_probabilities_adiabatic() takes the two tolerances as patch_atol and phase_tol, 1e-5 each by default (magnus.avgprob.PHASE_AVERAGE_PATCH_ATOL, magnus.avgprob.PHASE_AVERAGE_PHASE_TOL).

When there is no closed form

A profile with discontinuities — the PREM layer boundaries an Earth-crossing trajectory steps through — cannot be followed adiabatically across its jumps, so neither construction above applies. There, average=True propagates the probability at energies across a window and averages the results (magnus.avgprob.averaged_probabilities_numerically()).

This returns a different quantity from the other two paths. They return the phase average, which weights phases and samples no energy. This route returns the probability averaged over one particular window, a top-hat of 41 energies from a flavor state at the start, and the answer depends on the width of that window. The default width, magnus.avgprob.AVG_DEFAULT_ENERGY_SPREAD, is 10%, the order of a real detector’s energy resolution. Every call on this route emits magnus.oscprob.PhaseAveragingWarning naming the width, the number of samples and the standard error of the resulting mean, so the figure never depends silently on a constant the user did not choose. Users with a known resolution should pass their own, through average_spread, which on this route is the half-width of the window as a fraction of the energy, strictly between 0 and 1. The number of samples is average_n_samples, 41 by default (magnus.avgprob.AVG_DEFAULT_N_SAMPLES). The standard error of the mean falls as the inverse square root of that number, and each sample costs a full propagation, so a target error \(\delta\) needs about \(N(s/\delta)^2\) samples, with \(N\) the samples used and \(s\) the standard error the warning reported. The other routes sample nothing and refuse average_n_samples.

import numpy as np
import magnus.globaldefs as gd
import magnus.oscprob as oscprob

osc = gd.load_nufit_params('NuFIT 6.1')
edges = [1000.0*gd.UNIT_KM]                      # a density jump at 1000 km

def ne(l):                                       # electron density [eV^3]
    return np.where(l < edges[0], 100.0, 10.0)*gd.N_AV*gd.UNIT_PER_CM3

P = oscprob.osc_prob_matter_std_potential(
    3, ne, 10.0*gd.UNIT_MEV, 2000.0*gd.UNIT_KM, osc, L0=0.0, t_breakpoints=edges,
    average=True, average_spread=0.05, average_n_samples=161,
    density_is_of_number_of_electrons=True)

The density function may be written for one position or for an array of them, as ne is. Both give the same numbers; the array form is faster, since the engine samples the density at every quadrature node. Written for one position at a time, as 100.0 if l < edges[0] else 10.0, it issues a ScalarHamiltonianWarning that names rho_func.

Changed in version 1.2.0: average_spread sets the window on this route, and average_n_samples was added.

Cost

Case	Method	Cost
Vacuum, constant density (and their NSI/LIV variants)	Closed form	10 µs per energy where every phase has decohered, 23 µs where some has not
Exponential density, Sun (and their NSI/LIV variants), no window	Adiabatic transport	6 ms per energy
A profile with windows whose phases survive: a solar chord, 10 GeV to 1 TeV	Windows at Gauss–Hermite nodes	0.21 s (10 GeV) to 0.08 s (1 TeV) per point, at the default tolerance
Earth (PREM)	Sampled over an energy window	~0.1 s per energy

All rows were measured on an idle machine with NuFIT 6.1: the solar rows on 'BS05-AGS-OP', and the chords on 'B16-GS98' at \(b = 0.6\,R_\odot\). For comparison, averaging the engine’s output numerically over 2001 energies for the vacuum row takes about 0.25 s, and is itself an approximation.

Usage

import numpy as np

import magnus.oscprob as oscprob
import magnus.globaldefs as gd

# load_nufit_params returns just the six mixing parameters; the
# OSC_PARAMS_PREDEFINED entries also carry 'name'/'description' strings,
# which the propagation machinery would reject.
osc = gd.load_nufit_params('NuFIT 6.1')

# 1 TeV from 100 Mpc: every phase has decohered, and this is the limit
P = oscprob.osc_prob_3nu_vacuum(1.0*gd.UNIT_TEV, 3.0857e21*gd.UNIT_KM,
                                average=True, **osc)
np.round(np.asarray(P), 4)

array([[0.5481, 0.2155, 0.2364],
       [0.2155, 0.4071, 0.3774],
       [0.2364, 0.3774, 0.3862]])

The most quoted consequence of averaged astrophysical oscillations follows in one line: a source producing the pion-decay composition \((1:2:0)\) delivers something close to equipartition at Earth.

at_source = np.array([1.0, 2.0, 0.0])/3.0
at_earth = at_source @ np.asarray(P)
np.round(at_earth*3.0, 3)

array([0.979, 1.03 , 0.991])

At a terrestrial baseline, the phases have not decohered, and the answer depends on the spread:

import warnings

with warnings.catch_warnings():
    warnings.simplefilter('ignore')     # PhaseAveragingWarning: it depends on the spread
    P_mumu = [float(oscprob.osc_prob_3nu_vacuum(1.0*gd.UNIT_GEV, 1000.0*gd.UNIT_KM,
                    nu_i=gd.NUMU, nu_f=gd.NUMU, average=True, average_spread=s, **osc))
              for s in (0.05, 0.1, 0.3)]
np.round(P_mumu, 4)

array([0.9722, 0.9088, 0.5885])

The keyword is not confined to the wrappers. A Hamiltonian of your own goes through osc_prob_energy_baseline, or through osc_prob_earth and osc_prob_sun, which add the geometry. There, average=True takes the same three routes as on the wrappers. It uses the closed form when the Hamiltonian does not depend on position, and adiabatic transport when it does and the profile is smooth. It averages over an energy window when t_breakpoints or t_slab_edges declare discontinuities. The Hamiltonian is passed as H(E, l), or as H(E) with H_func_is_function_only_of_energy=True, or as a matrix. A matrix, or a function of position alone, does not depend on energy, so no spread acts on it; see above.

import magnus.hamiltonians as hamiltonians
import magnus.matter as matter

h_vac = hamiltonians.hamiltonian_3nu_vacuum_energy_independent(**osc)
proj = matter.matter_potential_projector(3)
vcc = matter.vcc_func_from_rho_func(
    lambda l: 5.0*np.exp(-l/(1000.0*gd.UNIT_KM)), 0.0, 1.0, 0.5,
    nubar=False, density_matter_is_in_g_per_cm3=True,
    density_is_of_number_of_electrons=False)

def H(E, l):
    return h_vac/E + np.asarray(vcc(l))[..., None, None]*proj

with warnings.catch_warnings():
    warnings.simplefilter('ignore')     # this one, too, depends on the spread
    P = oscprob.osc_prob_energy_baseline(
        H, 0.1*gd.UNIT_GEV, 5000.0*gd.UNIT_KM, 0.0,
        nu_i=gd.NUE, nu_f=gd.NUE, average=True)
round(float(P), 4)

0.1989

osc_prob itself, which computes one point, does not accept the keyword and raises ValueError if given it.

Checking whether a scan resolves the oscillation: `strategy_info['sampling']`

A scan of instantaneous probabilities over a long trajectory returns correct numbers that can still be the wrong quantity: if no detector resolves the phase, the observable is the average.

Every osc_prob_* entry point that accepts strategy_info reports how coarsely the request samples the oscillation it is computing:

R_sun = gd.SUN_RADIUS*gd.UNIT_KM
info = {}
P = oscprob.osc_prob_3nu_sun(10*gd.UNIT_MEV, np.linspace(0.1, 0.9, 5)*R_sun, 0.0,
                             nu_i=gd.NUE, nu_f=gd.NUE, strategy_info=info)
info['sampling']
# {'oscillation_length': 5.00e+10,   'cycles_over_trajectory': 6.34e+04,
#  'spacing': 7.05e+14,              'cycles_per_step': 1.41e+04,
#  'nyquist_points': 126808,         'aliased': True}

cycles_per_step is the number to read. Above about 0.5, the scan takes fewer than two samples per oscillation, so the returned array cannot represent the oscillation and must not be plotted or interpolated as a curve. The individual values are correct, but the curve through them is an artifact. nyquist_points says how many baselines would be needed to sample it properly.

On the physically motivated profile families in docs/dev/adversarial_batteries/, the counts are large:

trajectory	oscillations across it	baselines for Nyquist
Earth chord	~430	861
Solar, one scale height	~2200	4390
Supernova ray	~37 000	73 392

This is reported, not warned about. The Nyquist criterion would fire on 44 of 45 realistic scan sizes. A warning on 98% of calls would be noise, and would teach users to silence a category that also carries real discontinuity warnings. The measurement is adversarial_batteries/alias_fp.py.

The report requires eigenvalues at eight points along the trajectory, so it is computed only when strategy_info is supplied. Without it, the report costs nothing; with it, it costs 5.5% of the cheapest scan measured and under 0.1% of a substantial one.

When aliased is True, consider whether the average is the intended observable. If it is, average=True gives it without resolving the oscillation, with the spread set by average_spread, and magnus.oscprob.PhaseAveragingWarning says where the answer depends on that spread.

How much does the phase matter?

It depends on the profile. The table below compares the error of an instantaneous scan with that of its average over six oscillation lengths, both against a solve_ivp reference (adversarial_batteries/avg_check.py and avg_check2.py). The solar row uses the log-linear interpolant of the BS05 table, the profile of notebook 13; a cubic-spline interpolant gives 8.889e-4 and 6.051e-4, with the same conclusion:

configuration	instantaneous	averaged	averaged inside 1e-3
Solar model, 2 flavors, 5 MeV	6.000e-04	7.110e-04	yes
Supernova turbulence, 45 MeV	5.584e-03	3.843e-04	yes
Supernova shock, 70 km front	4.917e-04	2.151e-04	yes
Supernova shock, 0.07 km front	1.988e-01	2.222e-01	no, and warned

Only in the last row is the observable wrong. A shock front changes the adiabaticity of the level crossing, so it moves the conversion probability itself, not only the phase of its oscillation. That is an error in the amplitude, and no averaging removes it. In the other rows, the averaged answer is within the target even where a single baseline is not, because the instantaneous error is mostly in the phase: the profile shifts when the oscillation peaks, which no observable resolves.

Warning

Do not read the ratio of these two columns as a diagnostic. The averaged column is a finite-window mean, and such a mean is an estimator with a bias of its own. On a profile whose density varies across the averaging window, the bias does not shrink as the window widens, because a wider window also averages over different matter conditions. On the solar ray above, the window mean moves from 0.5962 to 0.6054 between six and forty-eight oscillation lengths, drifting away from a limit rather than toward one. Notebook 13 prints that sweep. The ratio is meaningful only on a controlled comparison at fixed matter conditions, as in notebook 23.

To obtain the averaged probability, ask for it rather than estimate it. average=True evaluates it with no window to choose. The solar ray has no non-adiabatic window, so there the phase average is the decohered limit, and it reproduces the adiabatic MSW expression \(\langle P_{ee}\rangle = \tfrac12 + \tfrac12\cos2\theta_m(L_0)\cos2\theta_m(L_1)\) to 3e-16 across 1–20 MeV.

Limitations and scope

The phase average is a statement about a measurement with an energy resolution. It is not a model of quantum decoherence: there is no dissipative term here. See the “When is Magνs not the right tool?” section of Magνs: neutrino oscillations via the Magnus expansion.
A spread in geometry — the size of the production region, the range of impact parameters a pixel covers — is not modeled. A phase that does not depend on energy, as the matter phases on a solar chord at TeV energies do, is kept whatever its size, and the user must average over the geometry separately.
The Earth/PREM path is a windowed average of the probability, as described above.
Both eigenbases, at production and at detection, come from the one Hamiltonian given. The limit of Position-dependent Hamiltonians also describes mixing that differs between the two ends for other reasons — a non-unitary mixing matrix, or parameters that run with energy — with \(P^{\rm cross}\) the identity and the two matrices supplied directly. Magνs does not accept them in that form.

See Available probability functions and the API reference for the full listing of magnus.avgprob.