magnus.avgprob

avgprob.py

Contains the phase-averaged (fully decohered) oscillation probabilities, the exact \(L/E \to \infty\) limit reached by high-energy astrophysical neutrinos.

Physical idea: a neutrino produced at a cosmological distance arrives with an oscillation phase \(\Delta m^2 L / 2E\) of order \(10^{15}\) or more, and neither the source distance, nor the production region, nor the detector’s energy resolution is known to anything close to that precision. Every oscillatory term is therefore averaged over many cycles and vanishes, leaving only the incoherent sum

\[P(\nu_\alpha \to \nu_\beta) = \sum_i |V_{\alpha i}|^2 |V_{\beta i}|^2 ,\]

where \(V\) diagonalizes the Hamiltonian. This is not an approximation to be refined: it is the exact limit, and it costs one matrix product rather than an integration. For standard vacuum oscillations the result does not depend on energy or baseline at all, so a single matrix serves an entire flux calculation.

Coherence is decided physically, not numerically

The formula above assumes every relative phase averages away. That is a statement about pairs of eigenvalues, not about the spectrum as a whole: the pair \((i,j)\) decoheres only if \((\lambda_i - \lambda_j) L\) sweeps through many cycles across the averaging window. Two eigenvalues that are close enough to keep their relative phase fixed stay coherent, and their cross term survives.

This module therefore groups the spectrum into blocks of mutually coherent eigenvalues and sums coherently inside each block,

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

which reduces to the familiar expression when every block is a singleton. The distinction is not academic here: a sterile state with a small \(\Delta m^2_{41}\), or any degenerate spectrum, makes the naive sum quietly wrong.

The same per-pair phase decides whether the averaged limit applies at all. A pair whose phase spread is neither much larger than \(2\pi\) (decohered) nor much smaller than one (coherent) sits in between, where no closed form is valid; coherence_report() names those pairs, and the callers in magnus.oscprob warn rather than return a number the physics does not support.

This module is self-contained: it depends only on numpy, not on magnus.oscprob, so it can be applied to any Hermitian Hamiltonian of any dimension independently of the rest of the API.

Routine listings

coherence_blocks coherence_report averaged_probabilities_from_eigenbasis averaged_probabilities_constant_hamiltonian adiabatic_phase_differences level_crossing_matrix averaged_probabilities_adiabatic averaged_probabilities_numerically

Attributes

DECOHERENCE_PHASE_THRESHOLD	Module-level constant
COHERENCE_PHASE_THRESHOLD	Module-level constant
AVG_DEFAULT_ENERGY_SPREAD	Module-level constant
AVG_DEFAULT_N_SAMPLES	Module-level constant

Functions

coherence_blocks(→ List[List[int]])	Groups eigenvalues into blocks that stay mutually coherent.
coherence_report(→ Tuple[List[List[int]], ...)	Reports the coherence structure of a spectrum, and which pairs sit in
averaged_probabilities_from_eigenbasis(→ numpy.ndarray)	Phase-averaged oscillation probabilities from the eigenbasis of the
averaged_probabilities_constant_hamiltonian(...)	Phase-averaged oscillation probabilities for a constant Hamiltonian.
averaged_probabilities_numerically(...)	Averages a probability by sampling it across an energy window.
adiabatic_phase_differences(→ numpy.ndarray)	Relative phases accumulated between instantaneous eigenvalues.
level_crossing_matrix(→ Tuple[numpy.ndarray, ...)	Probability of ending on level \(j\) having started on level \(i\).
averaged_probabilities_adiabatic(...)	Phase-averaged probabilities for a position-dependent Hamiltonian.

Module Contents

magnus.avgprob.DECOHERENCE_PHASE_THRESHOLD

Module-level constant

Accumulated phase spread, in radians, above which a pair of eigenvalues is treated as fully decohered. One full cycle is the point at which the average of \(\cos\Delta\phi\) over the window has collapsed to a small fraction of its coherent value, and every further cycle only reduces it.

Added in version 1.0.0.

Type:

float

magnus.avgprob.COHERENCE_PHASE_THRESHOLD = 0.01

Module-level constant

Accumulated phase spread, in radians, below which a pair of eigenvalues is treated as fully coherent, so that its cross term is kept in full.

The gap between this and DECOHERENCE_PHASE_THRESHOLD is deliberate and is not a tolerance to be tightened away: a pair falling between the two is in neither limit, and no averaged expression describes it. Such pairs are reported by coherence_report() rather than silently assigned to one side.

Added in version 1.0.0.

Type:

float

magnus.avgprob.coherence_blocks(eigenvalues: Sequence[float] | numpy.ndarray, phase_scale: float, decoherence_threshold: float | None = DECOHERENCE_PHASE_THRESHOLD) → List[List[int]]

Groups eigenvalues into blocks that stay mutually coherent.

Two eigenvalues belong to the same block when the phase they accumulate relative to each other, \(|\lambda_i - \lambda_j| \times\) phase_scale, stays below decoherence_threshold, so that their cross term in the probability is not averaged away.

Grouping is by transitive closure over that relation, which is the conservative choice: a chain of individually-close eigenvalues is kept in one block rather than split at an arbitrary point. A spectrum whose spacings are all comparable to the threshold therefore collapses into a single block, and is exactly the case coherence_report() flags as having no valid averaged limit.

Added in version 1.0.0.

Parameters:

• eigenvalues (list or np.ndarray) – Eigenvalues of the Hamiltonian [eV]. Need not be sorted.

• phase_scale (float) – Baseline over which the phase accumulates [\(\text{eV}^{-1}\)], so that (lambda_i - lambda_j)*phase_scale is a phase in radians.

• decoherence_threshold (float, optional) – Phase above which a pair is treated as decohered. Default: DECOHERENCE_PHASE_THRESHOLD.

Returns:

Indices of eigenvalues, grouped into blocks and sorted within each block. The blocks themselves are ordered by their smallest index, so the result is deterministic.

Return type:

list of list of int

Examples

A spectrum whose splittings are all large is fully decohered, one index per block; two eigenvalues sharing a value stay together.

import magnus.avgprob as ap

ap.coherence_blocks([0.0, 1.0, 2.0], phase_scale=1.0e3)

[[0], [1], [2]]

magnus.avgprob.coherence_report(eigenvalues: Sequence[float] | numpy.ndarray, phase_scale: float, decoherence_threshold: float | None = DECOHERENCE_PHASE_THRESHOLD, coherence_threshold: float | None = COHERENCE_PHASE_THRESHOLD) → Tuple[List[List[int]], List[Tuple[int, int, float]]]

Reports the coherence structure of a spectrum, and which pairs sit in neither limit.

Every pair of eigenvalues is in one of three regimes, set by the phase it accumulates relative to the others over phase_scale:

• far above decoherence_threshold, the cross term has averaged away and the pair contributes incoherently;

• far below coherence_threshold, the relative phase has barely advanced and the pair is still fully coherent;

• in between, neither statement holds, and no averaged expression is a valid description – the honest answer there is the full oscillation probability, not an average.

Added in version 1.0.0.

Parameters:

• eigenvalues (list or np.ndarray) – Eigenvalues of the Hamiltonian [eV].

• phase_scale (float) – Baseline over which the phase accumulates [\(\text{eV}^{-1}\)].

• decoherence_threshold (float, optional) – Phase above which a pair counts as decohered. Default: DECOHERENCE_PHASE_THRESHOLD.

• coherence_threshold (float, optional) – Phase below which a pair counts as fully coherent. Default: COHERENCE_PHASE_THRESHOLD.

Returns:

The coherence blocks, and the list of (i, j, phase) triples for pairs that are in neither limit. An empty second element means the averaged result is exact for this spectrum and baseline.

Return type:

(list of list of int, list of (int, int, float))

magnus.avgprob.averaged_probabilities_from_eigenbasis(eigenvectors: Sequence | numpy.ndarray, blocks: List[List[int]] | None = None) → numpy.ndarray

Phase-averaged oscillation probabilities from the eigenbasis of the Hamiltonian.

Computes

\[P_{\alpha\beta} = \sum_b \Big| \sum_{i \in b} V^*_{\alpha i} V_{\beta i} \Big|^2 ,\]

the sum over coherence blocks b of the squared modulus of the coherent amplitude within each block. With one index per block this is the familiar \(\sum_i |V_{\alpha i}|^2 |V_{\beta i}|^2\).

The result is symmetric, so the averaged probability is the same in both directions, and identical for neutrinos and antineutrinos: conjugating \(V\) leaves every term unchanged. CP violation does not survive the average, even though the mixing angles and phases do enter through \(|V_{\alpha i}|\).

Added in version 1.0.0.

Parameters:

• eigenvectors (list or np.ndarray) – Matrix whose columns are the eigenvectors of the Hamiltonian, shape (..., d, d). A leading batch axis is allowed and is broadcast over, so an array of energies costs one contraction.

• blocks (list of list of int, optional) – Coherence blocks, as returned by coherence_blocks(). If None (default), every eigenvalue is assumed to have decohered from every other, which is the astrophysical case.

Returns:

Averaged probability matrix, shape (..., d, d), with the initial flavor as the row index, so each row sums to one.

Return type:

np.ndarray

Examples

import numpy as np

import magnus.avgprob as ap
import magnus.hamiltonians as hams

U = hams.pmns_mixing_matrix(0.55, 0.68, 0.15, 3.7)
P = ap.averaged_probabilities_from_eigenbasis(U)
np.round(P, 4)

array([[0.5528, 0.2153, 0.2319],
       [0.2153, 0.4072, 0.3775],
       [0.2319, 0.3775, 0.3906]])

magnus.avgprob.averaged_probabilities_constant_hamiltonian(hamiltonian: Sequence | numpy.ndarray, baseline: float | None = None) → numpy.ndarray

Phase-averaged oscillation probabilities for a constant Hamiltonian.

Diagonalizes hamiltonian and applies averaged_probabilities_from_eigenbasis(). This covers every position-independent case – vacuum, matter of constant density, and their NSI and LIV variants – exactly, at the cost of one eigendecomposition.

Added in version 1.0.0.

Parameters:

• hamiltonian (list or np.ndarray) – Hermitian Hamiltonian [eV], shape (..., d, d). A leading batch axis (energies, say) is allowed.

• baseline (float, optional) – Baseline [\(\text{eV}^{-1}\)], used only to decide which eigenvalues have decohered from each other. If None (default), every pair is taken to be decohered, which is the astrophysical limit and makes the result independent of distance.

Returns:

Averaged probability matrix, shape (..., d, d), rows summing to one.

Return type:

np.ndarray

magnus.avgprob.AVG_DEFAULT_ENERGY_SPREAD = 0.1

Module-level constant

Half-width of the energy window, as a fraction of the energy, used when the averaged probability has to be obtained by sampling rather than in closed form.

Ten per cent is the order of a real detector’s energy resolution, and it is the smearing that does the averaging: the physical statement is that the oscillation phase varies by many cycles across whatever window the measurement integrates over. It is a default, not a property of the physics, so it is named here rather than buried, every use of it is warned about, and callers with an actual resolution should pass theirs.

Added in version 1.0.0.

Type:

float

magnus.avgprob.AVG_DEFAULT_N_SAMPLES = 41

Module-level constant

Number of samples across the window used by averaged_probabilities_numerically().

The sampled phases are effectively independent when the accumulated phase is large, so the error of the mean falls only as \(1/\sqrt{N}\) – 41 samples give a few per cent. Raising it buys accuracy slowly and costs a full propagation each; the closed-form paths in this module exist precisely to avoid this trade.

Added in version 1.0.0.

Type:

int

magnus.avgprob.averaged_probabilities_numerically(prob_of_energy: Callable, energy: float, relative_spread: float | None = AVG_DEFAULT_ENERGY_SPREAD, n_samples: int | None = AVG_DEFAULT_N_SAMPLES) → Tuple[numpy.ndarray, float]

Averages a probability by sampling it across an energy window.

The fallback for cases with no closed form – a profile with discontinuities, say, where there is no instantaneous eigenbasis to decohere in. Unlike the closed forms in this module, this is not the \(L/E \to \infty\) limit: it is the average over a particular window, and the answer depends on that window. Its width is therefore an argument, and callers that leave it at the default should say so to their own callers.

Samples are uniform in \(1/E\), in which the oscillation phase is linear, so they are spread evenly in phase rather than bunched.

Added in version 1.0.0.

Parameters:

• prob_of_energy (Callable) – Returns the probability matrix at a given energy; called once per sample.

• energy (float) – Central energy [eV].

• relative_spread (float, optional) – Half-width of the window as a fraction of energy. Default: AVG_DEFAULT_ENERGY_SPREAD.

• n_samples (int, optional) – Number of samples. Default: AVG_DEFAULT_N_SAMPLES.

Returns:

The mean probability matrix, and the largest standard error of the mean across its entries – the honest uncertainty of the result, which a closed form would not have.

Return type:

(np.ndarray, float)

magnus.avgprob.adiabatic_phase_differences(H_func: Callable, l0: float, l1: float, n_points: int | None = 201) → numpy.ndarray

Relative phases accumulated between instantaneous eigenvalues.

In the adiabatic regime a neutrino stays on one level and accumulates the dynamical phase \(\int \lambda_i(l)\, dl\), so the phase that decides whether levels \(i\) and \(j\) still interfere is \(\Delta\phi_{ij} = \int_{l_0}^{l_1} [\lambda_i(l) - \lambda_j(l)]\, dl\). That integral, not the eigenvalue gap at any single point, is what the coherence tests in this module need for a position-dependent Hamiltonian.

Integrated with Simpson’s rule: the trapezoid leaves a residual here that is easily mistaken for a physical effect (the same error, in the same integral, once looked like a floor on the accuracy of adiabatic transport in magnus.adiabatic).

Added in version 1.0.0.

Parameters:

• H_func (Callable) – Hamiltonian as a function of position, H_func(l) [eV].

• l0 (float) – Start and end of the trajectory [\(\text{eV}^{-1}\)].

• l1 (float) – Start and end of the trajectory [\(\text{eV}^{-1}\)].

• n_points (int, optional) – Number of sampling points; forced to be odd for Simpson’s rule. Default: 201.

Returns:

Matrix of accumulated phase differences, shape (d, d), antisymmetric.

Return type:

np.ndarray

magnus.avgprob.level_crossing_matrix(H_func: Callable, l0: float, l1: float, threshold: float | None = 0.1, n_probe: int | None = 200, fd_step_frac: float | None = 1e-06, magnus_exp_order: int | None = 6, integration_method: str | None = 'gl') → Tuple[numpy.ndarray, List[Tuple[float, float]], bool]

Probability of ending on level \(j\) having started on level \(i\).

Adiabatic evolution keeps a neutrino on the level it was produced on, so this matrix is the identity wherever the adiabatic approximation holds. It departs from the identity only across a non-adiabatic window – a resonance sharp enough for levels to exchange character faster than the state can follow – and it is exactly there that the averaged probability needs it.

The window is located with the Hellmann-Feynman diagnostic in magnus.adiabatic, and the transfer across it is computed with that module’s own convergence-checked Magnus patch rather than with a Landau-Zener formula, so it inherits an exact treatment of the crossing instead of an asymptotic approximation to it.

Added in version 1.0.0.

Parameters:

• H_func (Callable) – Hamiltonian as a function of position, H_func(l) [eV].

• l0 (float) – Start and end of the trajectory [\(\text{eV}^{-1}\)].

• l1 (float) – Start and end of the trajectory [\(\text{eV}^{-1}\)].

• threshold (float, optional) – Adiabaticity threshold passed to magnus.adiabatic.find_nonadiabatic_windows(). Default: 0.1.

• n_probe (int, optional) – Density of the search grid for the same. Default: 200.

• fd_step_frac (float, optional) – Finite-difference step, as a fraction of the domain, for the same. Default: 1e-6.

• magnus_exp_order (int, optional) – Magnus order for the local patch. Default: 6.

• integration_method (str, optional) – Integration method for the local patch. Default: ‘gl’.

Returns:

The level-to-level probability matrix, with the starting level as the row index; the non-adiabatic windows found; and whether every local patch converged. A False in the last position means the crossing probabilities are not trustworthy, not that they are merely imprecise.

Return type:

(np.ndarray, list of (float, float), bool)

magnus.avgprob.averaged_probabilities_adiabatic(H_func: Callable, l0: float, l1: float, n_points: int | None = 201, threshold: float | None = 0.1, n_probe: int | None = 200, fd_step_frac: float | None = 1e-06, magnus_exp_order: int | None = 6, integration_method: str | None = 'gl') → Tuple[numpy.ndarray, dict]

Phase-averaged probabilities for a position-dependent Hamiltonian.

A neutrino produced at \(l_0\) decoheres in the eigenbasis there, is carried along the levels of the instantaneous Hamiltonian, and is detected in the eigenbasis at \(l_1\):

\[P_{\alpha\beta} = \sum_{ij} |V_{\alpha i}(l_0)|^2\, P^\text{cross}_{ij}\, |V_{\beta j}(l_1)|^2 ,\]

with \(P^\text{cross}\) from level_crossing_matrix() – the identity wherever the evolution is adiabatic. This is the standard MSW-plus-decoherence result, generalized to any number of levels and any number of crossings.

Two things have to hold for the expression to mean anything, and both are checked rather than assumed. The levels must have decohered from each other by the time of detection, and if there is more than one crossing they must also have decohered between crossings, since otherwise composing the crossings as probabilities – rather than as amplitudes – discards interference that is still there. Both are reported.

Added in version 1.0.0.

Parameters:

• H_func (Callable) – Hamiltonian as a function of position, H_func(l) [eV].

• l0 (float) – Production and detection positions [\(\text{eV}^{-1}\)].

• l1 (float) – Production and detection positions [\(\text{eV}^{-1}\)].

• n_points (int, optional) – Sampling density for the accumulated-phase integrals. Default: 201.

• threshold (float, int, float, optional) – Passed to level_crossing_matrix().

• n_probe (float, int, float, optional) – Passed to level_crossing_matrix().

• fd_step_frac (float, int, float, optional) – Passed to level_crossing_matrix().

• magnus_exp_order (int, optional) – Magnus order for the local patches. Default: 6.

• integration_method (str, optional) – Integration method for the local patches. Default: ‘gl’.

Returns:

The averaged probability matrix, rows summing to one, and a report with keys 'windows' (the non-adiabatic windows), 'patches_converged' (bool), 'undecided' (pairs that are in neither the coherent nor the decohered limit over the whole trajectory) and 'undecided_between_crossings' (the same, over each adiabatic stretch separating two crossings).

Return type:

(np.ndarray, dict)